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1.  INTRODUCTION 


Unsteady  flow  in  the  passages  of  wave  rotor  devices  can  adequately  be 
modelled  on  a  one -dimensional  basis.  However,  this  modelling  can  be  quite 
involved  due  to  the  peculiar  characteristics  typical  of  wave  rotor  type  flows. 
The  numerical  calculation  has  to  provide  approximate  solutions  of 
time-dependent  compressible  fluid  flow  problems  which  involve  discontinuities 
and  strong  wave  interactions.  Ref.  (1)  lists  three  criteria  which  such 
approximate  solutions  should  satisfy  simultaneously:  (i)  the  solution  must  be 
reasonably  accurate  in  smooth  regions  of  the  flow.  Continuous  waves 
(rarefaction  waves,  compression  waves)  should  propagate  at  the  correct  speed 
and  should  maintain  the  correct  shape  which  involves  steepening  or  spreading 
at  the  correct  rate;  (ii)  discontinuities  which  are  transported  along 
characteristics  (gradient  discontinuities,  contact  surfaces),  should  be 
modelled  by  sharp  and  discrete  jumps,  and  should  be  transported  at  the  correct 
speed;  and  (iii)  nonlinear  discontinuities  such  as  shocks  should  be  computed 
stably  and  accurately. 

In  addition,  the  complex  pattern  of  shock  waves  and  contact  surfaces  that 
could  evolve  in  wave  rotor  devices  precludes  the  use  of  numerical  methods 
which  rely  on  either  some  type  of  artificial  viscosity  or  a  special  treatment 
of  discontinuities.  Such  methods  would  quickly  become  quite  impractical  for 
this  application  due  to  programming  difficulties  and  cost  of  execution. 

Computation  of  such  solutions  has  generally  been  carried  out  by  solving  a 
set  of  finite  difference  equations  which  approximate  the  governing 
differential  equations  of  flow.  All  such  schemes  inherently  have  a  finite 
amount  of  dissipation  as  well  as  dispersion  of  the  wave  modes  they  model,  and 
it  is  difficult  to  construct  difference  schemes  which  simultaneously  satisfy 
the  criteria  given  above.  Stability  problems  may  also  be  an  added  concern  for 


these  schemes 


In  view  of  the  foregoing,  an  alternative  approach  to  solving  wave  rotor 
type  flows  was  sought,  and  the  purpose  of  this  report  is  to  describe  such  a 
scheme  along  with  some  results.  The  scheme  is  known  variously  as  Glimm's 
method,  the  Random  Choice  Method  (RCM)  or  the  piecewise  sampling  method.  The 
method  evolved  from  a  constructive  proof  of  the  existence  of  solutions  to 
systems  of  nonlinear  hyperbolic  conservation  laws  given  by  Glimm  (Ref.  2). 
Chorin  (Refs.  3  and  4)  developed  the  scheme  into  an  effective  numerical  tool 
for  gas  dynamic  applications,  with  emphasis  on  detonation  combustion  problems 
and  reacting  gas  flows.  Although  the  RCM  computes  solutions  on  a  fixed  grid, 
it  is  not  a  difference  scheme,  utilizing  solutions  of  locally  defined  Rlemann 
problems  as  the  basic  building  blocks  for  the  global  solution.  Each  of  the 
local  Rlemann  problems  (defined  in  more  detail  in  section  2)  provides  an 
analytically  exact  elementary  similarity  solution.  By  means  of  a  suitable 
sampling  procedure,  usually  of  a  pseudo-random  or  quasi-random  nature,  the 
similarity  solutions  are  superposed  to  construct  the  approximate  solution  to 
the  equations. 

With  an  appropriate  sampling  technique,  the  RCM  in  one  dimension  is 
possibly  superior  to  any  finite  difference  scheme  in  meeting  the  criteria 


established  above 


2.  METHOD 


2.1.  Solution  Procedure 

The  method  models  the  one -dimensional,  compressible,  inviscid  Euler 

equations,  expressed  in  conservation  form  as 

iU  ^  3F(U).  .  0  ,  where 

3t  3x 


U(x,t)  *  jpu 


F(U)  =  pu2  +  p 
(E  +  p)p 


Here  E  is  the  total  energy  per  unit  volume  and  may  be  expressed  as  (for  a 
poly tropic  gas) 

E  »  ,je  +  —  pu2  ,  e  J  internal  energy  per  unit  mass 

■7^  (a) 

p  is  the  density,  p  is  pressure  and  u  is  velocity  in  the  one  space 
dimension  being  considered  here.  With  initial  data  specified  in  the  form 

U(x,0)  =  9(x)  , 

an  initial  value  problem  is  defined  for  the  Euler  equations.  The  simplest 
initial  value  problem  for  which  discontinuities  appear  is  the  Riemann  problem 
to  find  the  gas  flow  resulting  from  an  initial  state  in  which  the  gas  on  the 
right  of  an  'origin'  is  in  a  constant  state,  and  the  gas  on  the  left  is  in 


another  constant  state,  i.e. 


X  <  0 


X  >  0 


V'L,R 

.,R  “  (►Ju)l,R 

El,r 


3 


where  the  subsripts  L  and  R  denote  the  left  and  right  sides  of  the 


'origin',  here  arbitrarily  prescibed  at  0  .  That  is,  the  Riemann  problem 

consists  of  prescribing  constant  initial  data  on  either  side  of  an  origin 
where  a  jump  discontinuity  exists.  As  mentioned  before,  the  solution  of  the 
problem  constitutes  a  basic  building  block  of  the  random  choice  method.  A 
special  case  of  the  Riemann  problem  in  which  ul  =  ur  =  0  is  often  referred 
to  as  the  shock  tube  problem.  The  answer  to  the  problem  is  that  there  are 
four  possible  types  of  subsequent  flow,  depending  on  the  inequalities  in  the 
left  and  right  side  data  prescribed.  Thus,  in  both  directions  from  the 
origin,  a  shock  or  a  centered  rarefaction  wave  may  propagate,  giving  rise  to 
the  above  mentioned  four  different  possibilities.  Fig.  (1)  illustrates  the 
special  case  of  shock  tube  type  flow  and  the  evolution  of  the  wave  pattern. 

Fig.  (2)  shows  the  simple  fixed  Cartesian  grid  set  up  for  the  method. 

Let  Ax  be  a  spatial  increment  and  At  a  time  increment.  The  solution  is  to 
be  evaluated  at  time  (n  +  l)At  ,  n  being  a  non-negative  integer,  at 
spatial  increments  iAx  ,  i=  1,2, 3,  ...  The  initial  data  is  prescribed 
for  each  time  step  at  nAt  in  a  piecewise  constant  manner  i.e.,  it  consists 
of  intervals  of  length  Ax  where  the  data  is  constant,  separated  by  jump 
discontinuities : 


U(x  ,  nAt)  =  U?  ,  (i-^)Ax  <  x  <  (i+j)Ax 

The  solution  at  time  (n+l)At  then  is  required  to  have  the  same  property, 
i.e.,  it  is  piecewise  constant  over  an  interval  Ax  ,  and  it  serves  as  the 


initial  data  for  the  next  time  step: 

U(x,(n+l)At)  =  ,  (i  -  y)Ax  <  X  <  (i+y)Ax 


This  procedure  defines  a  sequence  of  local  Riemann  problems  to  be  solved  at 
each  time  level.  On  the  grid  shown  in  Fig.  2,  for  example,  initial  data  wou 
be  specified  at  points  1,  3,  5  ....,  setting  up  a  succession  of  Riemann 
problems  defined  by  each  pair  of  states  (1,3)  ,  (3,5)  ,  (5,7),  with  the 
discontinuities  at  the  midpoint  of  each,  i.e.,  at  2,  4,  6,  ....  etc.  If  the 
time  step  increment  At  is  calculated  such  that 

At  <  0. (Ax),  max  (ju”|+a")  ,  with 

0  <  o  <  ^ 

then  the  waves  generated  at  the  discontinuities  of  adjacent  Riemann  problems 

will  not  interact,  as  shown  schematically  in  Fig.  2. 

Each  of  the  local  Riemann  problems  yields  an  exact  analytical  solution 

with  the  resulting  wave  structure  a  particular  combination/variation  of  the 

general  structure  shown  in  Fig.  3. 

In  the  x-t  plane,  the  solution  to  a  Riemann  problem  consists  of 

essentially  four  regions  connected  by  three  waves.  Thus  states  1  and  IV  are 

the  prescribed  left  and  right  states  for  the  problem,  and  states  II  and  III 

are  the  ’starred'  middle  states  separated  by  a  slip  line  or  contact 
dx 

discontinuity  =  u*  .  The  velocity,  u  ,  and  pressure,  p  ,  are 
continuous  across  the  contact,  but  tj  in  general  is  not.  Thus  ul*  =  ur* 

PL*  “  PR*  and  pl*  *  pr*  .  Si,b  ,  S2,b  and  Si,f  ,  S2,f  represent 

respectively  the  backward  and  forward  facing  waves  generated  at  the  point  of 

discontinuity  and  may  be  either  shocks  or  rarefaction  waves. 

Still  referring  to  Fig.  3,  it  is  seen  that  at  a  time  nAt  <  t  <  (n+l)At 

the  exact  solution  of  the  local  Riemann  problem  for  the  interval  [(i-l)Ax 

iAx]  may  actually  consist  of  several  distinct  states.  Consider  now  a 


translation  of  each  interval  [(i-l)^ix,  iAx]  to  ,  +  ±1^  J  such  that 

the  discontinuity  (i.e.,  the  point  from  which  the  waves  are  generated)  is 

centered  at  a  zero  origin.  Let  G  be  the  value  of  a  random  variable,  defined 
over  the  interval  [~-  »  +  >  ^nd  let 

C  =  OAx  ,  i.e. - 2  ^  ^  ^  ~2 

Also,  define  (x,t)  ,  nAt  <  t  <  (n+l)At  ,  to  be  the  exact  solution  to 

exact 

each  Riemann  problem.  Using  the  value  of  C  to  fix  a  point  in  the  interval 
Ax  of  each  Riemann  problem,  the  exact  solution  at  that  point  is  determined 
and  assigned  to  either  the  left  or  the  right  grid  point,  depending  on  whether 
^  is  <  or  >  0  .  Thus,  if  the  point  fixed  by  C  is  P’  (Fig.  3),  the 
exact  solution  to  the  Riemann  problem  at  that  sampled  location  is  assigned  to 
the  grid  point  on  the  right  and  if  the  sampled  point  is  P"  ,  the  solution  at 
that  location  is  assigned  to  the  grid  point  on  the  left,  i.e.,  for  a  typical 


interval  ((i-l)Ax,  iAx]  , 


5  <  0  .  =  U 


exact 


(C  ,  t) 


and  If  ?  >  0  .  uj-"*  .  (C  ,  t) 

It  is  seen  immediately  that  although  the  solutions  are  computed  on  a 
grid  in  this  method,  it  is  not  a  differencing  scheme.  Also,  instead  of  using 
a  weighted  average  of  the  Riemann  problem  solution  to  arrive  at  the  solution 
for  a  grid  pointt,  the  RCM  samples  a  particular  value  from  an  explicit  wave 
t  The  Godunov  method,  for  example  Implements 


n+l  1 

j  (1-T7)ux 


solution,  thus  eliminating  the  smoothing  out  of  wave  transport  and  Interaction 
Information  Inherent  In  averaging.  This  leads  to  the  'infinite*  resolution  of 
contact  discontinuities  and  shocks  that  the  scheme  displays. 

From  the  foregoing  discussion,  it  is  evident  that  the  success  of  the 
scheme  hinges,  to  a  large  extent,  on  the  inexpensive  and  exact  solution  of 
Riemann  problems  and  an  appropriate  sampling  technique.  Ref.  (3)  describes  a 
modification  to  an  Iterative  method  due  to  Godunov  (Ref.  5).  Theoretical 
details  for  the  Riemann  problem  solution  are  also  given  in  Ref.  (6). 

The  mathematical  properties  required  in  a  sampling  procedure  applicable 
to  this  scheme  are  defined  in  Ref.  (1).  A  brief  description  of  the  procedure 
is  given  below. 

In  previous  computations  using  the  RCM,  random  sampling  with  some 
variance  reduction  technique  (stratified  sampling),  was  used,  i.e.,  the  values 
were  taken  from  the  random  number  generator  installed  in  the  computer  (Ref. 

3).  It  was  shown  in  Ref.  (1)  that  a  more  accurate  form  of  sampling  is  a 
technique  due  to  van  der  Corput  (Ref.  7).  The  sequence  generated  is,  strictly 
speaking,  non-random,  but  has  particular  statistical  properties  that  are 
suitable  to  the  application.  The  sequence  is  referred  to  as  quasirandom  and 
is  generated  as  follows: 

The  binary  expansion  of  natural  numbers  may  be  expressed  as  (with  R=2): 

n  =  AqR®  +  AirI  +  A2r2  + . +  A^jR™  ,  (0  <  A  ^  <  R) 

m 

i.e.  n  =  I  Ak.z'^  ,  with  Aj^  =  0  or  1  ,  n  =  1,  2,  3,  ... 

k=0 


7 


Next,  the  digits  of  the  binary  numbers  are  reversed  and  a  decimal  point  is  put 
preceding  the  number;  this  gives  the  numbers 


<;>n  *  AqR"^  +  AiR“2  +  ...  + 


m 

9n  *  I  ,  again  with  Afc  ■  0  or  1 

k-0 


Conversion  to  the  decimal  scale  of  these  numbers  yields  the  required  sequence 

of  quasirandom  numbers  defined  over  the  interval  [0,1],  i.e., 

(t>n  (decimal)  ■  Gn  +  y 

®n  *  <i>n  (decimal)  - 
and  Cn  ••  On*Ax  as  defined  earlier. 


The  first  few  elements  of  the  sequence  given  below  illustrate  the  construction 


n*l  (decimal) 

*  1  (binary); 

<Pl  “  0.1  (binary) 

■  0.5  (decimal) 

2 

10 

0.01 

0.25 

3 

11 

0.11 

0.75 

4 

100 

0.001 

0.125 

5 

101 

0.101 

0.625 

6 

110 

0.011 

0.375 

7 

111 

0.111 

0.875 

8 

1000 

0.0001 

0.0625 

The  van  der  Corput  sequence  is  'equidistributed ' ,  and  yields  better  results 
than  those  obtained  using  a  'stratified'  random  sampling  technique. 

The  subroutine  employed  in  the  program  to  compute  the  random  numbers  is 
described  in  Appendix  B. 


2.2.  Boundary  Conditions 

In  general,  the  Implementation  of  boundary  conditions  in  the  RCM  is 
quite  straightforward,  but  does  require  some  thought.  Referring  to  Fig.  2, 
the  b.c.'s  are  specified  at  points  1  and  N  for  the  left  and  right  boundary 


respectively.  Note  that  if  the  sampled  solution  at  (n+l)iit  corresponds  to  a 
random  number  <  0  ,  the  solution  is  assigned  to  the  grid  point  on  the 
left.  For  the  Riemann  problem  defined  by  points  1  and  3,  the  sampled  solution 
would  then  be  assigned  to  grid  point  1  at  (n+l)At  ;  however  this  is 
overridden  by  assigning  the  proper  boundary  condition  at  1  again,  and  there  is 
no  contradiction.  A  similar  procedure  is  adopted  at  the  right  hand  boundary 
when  Cn  >  0  . 

The  subroutines  for  the  boundary  conditions  are  named  in  the  format 
BCxn  ,  BC  standing  for  boundary  Condition,  x  being  either  L  (for  ^eft), 
or  R  (for  Mght  boundary)  and  n  being  a  number  from  1  to  5  with  the 
following  designations: 

1  -  solid  wall  condition 

2  •*  outflow  at  constant  static  pressure 

3  -  special  formulation  ('piston'  inflow) 

4  -  isentropic  inflow  from  reservoir 

5  -  special  formulation  (rarefaction  wave  cancellation) 

2.2.1.  Solid  Wall  Conditions 

The  solid  wall  boundary  condition  requires  a  zero  normal  velocity  at 
the  wall  for  inviscid  flow  computations.  Due  to  the  random  sampling  involved 

A  V 

in  the  method  and  the  lateral  movement  of  the  sampled  solution  — ^  to  the 
left  or  right  of  the  discontinuity,  the  condition  is  difficult  to  Implement 

uniquely.  However,  the  procedure  adopted  here  is  found  to  yield  reasonably 

accurate  results  for  the  applications  intended.  (Note  that  the  difficulty  is 

not  unique  to  this  method  only.  The  implementation  of  zero  mass  flux  '.hrough 

a  surface  is  difficult  per  se  for  the  Euler  equations). 

Referring  to  Fig.  2,  let  the  physical  boundaries  be  at  point  2  and 


point  (N-1)  for  the  left  and  right  sides  respectively.  However,  the  boundary 
conditions  are  specified  at  point  1  (point  N)  for  the  left  (right)  side  as  a 
fictitious  'mirror'  state  of  the  conditions  at  point  3  (point  (n-2)) 
respectively,  but  with  the  reverse  sign  taken  for  the  velocity  component. 

Thus,  for  the  left  hand  boundary  Riemann  problem, 

“  P(3)  »  PL  “  P(3)  ,  UL  ■  -u(3) 

PR  =  P(3)  ,  PR  -  P(3)  ,  UR  *  u(3) 
and,  analogously,  for  the  right  hand  boundary  Riemann  problem. 

PL  *  p(N-2)  ,  PL  *  p(N-2)  ,  UL  *  u(N-2) 

PR  “  p(N-2)  ,  PR  =  p(N-2)  ,  UR  =  -u(N-2) 

The  solutions  are  then  sampled  in  the  manner  outlined  earlier. 

2.2.2.  Outflow  Conditions 

For  subsonic  outflow,  only  the  static  pressure  p  is  defined,  with 
the  continuation  condition  being  applied  to  the  rest  of  the  variables.  Thus, 
for  the  right  hand  boundary  for  example,  the  Riemann  problem  is  defined  as 
follows: 

PL  “  p(N-2)  ,  PL  =  p(N-2)  ,  UL  “  u(N-2) 

PR  “  Pout  •  PR  *  p(N-2)  ,  UR  =  -u(N-2) 
where  Pout  is  the  specified  outlet  pressure.  If  the  flow  going  out  is 
supersonic,  there  can  be  no  propagation  of  disturbances  upstream,  and  the 
continuation  condition  is  implemented  for  all  the  variables,  i.e.,  the  Riemann 
problem  now  is  the  trivial  case  defined  by 

PL  ”  p(N-2)  ,  PL  *  p(N-2)  ,  UL  *  u(N-2) 

Pr  *  p(N-2)  ,  PR  »  p(N-2)  ,  UR  «  -u(N-2) 

2.2.3.  Special  Formulation  of  'Piston'  Inflow 


In  general,  for  idealized  wave  rotor  flows,  hot  combustion  gases  are 


Introduced  into  the  rotor  through  nozzles  angled  such  as  to  allow  the  flow  to 
'slip  onto'  the  rotor,  i.e.,  without  incurring  incidence  or  deviation  angle 
losses.  Also,  in  the  ideal  treatment,  the  air  in  the  passages  of  a  wave  rotor 
is  exposed  to  the  hot  gas  at  high  pressure  instantaneously.  The  idealizations 
allow  for  uniform  conditions  to  be  prescribed  at  the  hot  gas  inlet  port. 

Thus,  a  'special'  form  of  inflow  boundary  condition  needs  to  be  specified 
here,  namely,  the  static  pressure,  the  velocity  and  the  density  of  the 
incoming  hot  gas.  Although  equivalent  to  specifying  the  total  pressure  and 
temperature  in  the  usual  Inflow  boundary  condition  treatment,  some  thought  is 
required  in  wave  rotor  type  flows  when  specifying  Pgas  *  ^^gas  ^nd  Ug^g 
This  is  because  only  a  shock  wave  needs  to  be  generated,  with  no  waves 
travelling  opposite  to  the  direction  of  flow.  The  solution  to  the  Riemann 
problem  would  then  consist  of  just  two  states  connected  by  a  single  shock 
wave.  The  flow  is  equivalent  to  that  generated  when  a  piston  is  pushed 
instantaneously  into  a  gas  at  rest.  In  general,  the  state  of  the  air  inside 
the  rotor  passage  is  known;  explicit  relations  for  two  states  connected 

through  a  shock  wave  are  given  in  Ref.  (6).  These  so-called  transition 

functions  help  in  specifying  the  boundary  conditions  for  the  incoming  flow 
properly. 

If  we  consider  the  left  boundary  for  this  inflow,  the  Riemann  problem 
is  set  up  as: 

PL  Phot  gas  •  ^^L  *  Phot  gas  t  ul  =  u^^ot  gas 

PR  “  P(3)  ,  PR  »  P(3)  ,  UR  “  u(3) 

with  pl  ,  PL  '*L  having  been  chosen  in  accordance  with  the 

considerations  discussed  above. 


2.2.4.  Isentropic  Inflow  From  Reservoir 


The  induction  of  fresh  charge  or  air  onto  the  rotor  usually 
corresponds  to  an  Isentropic  Inflow  situation.  The  flow  in  the  vicinity  of 
the  passage  end  can  be  treated  as  quasi-steady,  with  the  assumption  that  no 
flow  separation  takes  place  when  the  flow  enters.  Two  boundary  conditions  are 
required  for  this  type  of  Inflow;  these  are  provided  by  the  conservation  of 
energy  in  the  flow  from  the  external  region  to  the  inlet  (assumed  to  be 
steady),  and  by  the  presclbed  entropy  level  of  the  gas  in  the  external 
region. 

The  boundary  conditions  may  thus  be  expressed  as 

2^  2  2  2  2 

'^in  Y“1  ^in  Y~1  ^tot 

Sin  •  Stot 

where  the  subscripts  'in’  and  'tot'  apply  to  conditions  at  the  inlet  of  the 
passage  and  external  reservoir  respectively.  The  sonic  velocity  is  denoted  by 
a  ,  and  flow  velocity  by  u  .  Note  that  knowledge  of  the  Riemann  variable 
arriving  at  the  passage  end  from  within  the  passage  is  required  to  be  able  to 
solve  the  energy  equation  above  for  a^n  and  uin  .  For  the  left  end,  for 
example , 

„  2 

Qin  “  yTY  ®in  ~  '*in 

which  together  with  the  energy  equations  yields 


The  simple  analytical  treatment  given  above  has  to  be  modified 
somewhat  if  a  contact  discontinuity  is  formed  when  the  inflow  begins.  This 
is  due  to  the  fact  that  the  value  of  the  arriving  Riemann  variable  is  changed 
across  such  a  discontinuity,  which  thus  leads  to  an  additional  unknown. 
Procedures  for  solving  the  Inflow  for  these  situations  are  given  in  Ref.  (8). 
In  the  program  developed  here,  reasonably  good  results  are  obtained  by  setting 
the  velocity  at  the  boundary  point  equal  to  the  velocity  at  the  point  nearest 
the  physical  boundary.  For  the  left  end  e.g.,  the  variables  for  the  left 
state  of  the  Riemann  problem  are  obtained  as  follows: 

u(l)  “  u(3)  , 

a  reasonably  accurate  assumption  just  at  the  point  of  inlet  opening. 

Then,  from  the  'energy  ellipse'. 
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which  occupy  large  zones  in  the  passages  of  wave  rotors.  Fig.  (4)  shows  a 
wave  diagram  proposed  by  Spectra  Technology,  Inc.,  which  incorporates 
so-called  'wave  management'  or  'tuning'  ports  to  ideally  cancel  (and  otherwise 
attenuate)  impinging  rarefaction  fans.  The  physical  boundary  conditions  are 
thus  dictated  by  the  flow  developing  in  the  passage,  l.e.,  the  port  has 
non-uniform  flow  conditions  in  it,  which  at  each  point  match  those  of  the  flow 
at  the  end  of  the  passage  so  as  to  disallow  any  reflections  to  take  place. 
Numerically,  this  condition  is  achieved  by  implementing  the  continuity 
condition  across  the  boundary  for  all  the  flow  variables  involved.  For  the 
left  boundary,  thus,  the  Riemann  problem  is  defined  by: 


PL  ”  P(3) 

,  PL  *  P(3) 

,  UL  *  u(3) 

PR  *  P(3) 

,  PR  *  P(3) 

,  UR  -  u(3) 

and  analogously  for  the  right  boundary.  Note  that  these  boundary  conditions 
may  involve  either  inflow  or  outflow. 

2.3  Example  Calculations 

The  listing  of  the  program  is  included  in  Appendix  A,  and  the  various 
I  names  for  the  variables  are  listed  in  Appendix  B,  along  with  some  instructions 

on  how  to  use  the  program.  No  effort  as  yet  has  been  made  to  optimize  the 
code  either  for  storage  requirements  or  for  execution  efficiency. 

I 

I  In  this  section,  some  sample  calculations  are  carried  out  using  the  code, 

to  illustrate  its  usefulness  in  constructing  idealized  design  point  wave 
diagrams  which  can  serve  as  the  starting  configuration  for  detailed 
I  construction  of  diagrams  incorporating  real  flow  effects. 

2.3.1.  Test  Case  for  1-D,  Invlscld,  Unsteady,  Compressible  Flow 

fig*  (1)  illustrates  the  Initial  conditions  in  a  shock  tube,  with  the 

! 

'  diaphragm  at  xq  .  Sod  (Ref.  9)  suggested  a  test  case  for  hyperbolic 
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conservation  laws  with  the  following  conditions  as  initial  states  in  the  shock 
tube : 

PI  “  1.0  ,  PI  =  1.0  ,  ui  =  0.0 

P5  *  0.1  ,  P5  »  0.125  ,  U5  *  0.0 

i.e.,  the  gas  on  either  side  of  the  diaphragm  is  in  a  quiescent  state 
initially.  The  ratio  of  specific  heats  is  chosen  to  be  7/5,  and  Ax  is 
chosen  to  be  0.01. 

The  solution  (before  any  wave  has  reached  either  the  left  or  right 
end)  is  shown  in  Fig.  (5).  The  squares  shown  at  locations  xj  ,  X2  ,  X3 
and  X4  in  the  density  plot  give  the  analytically  calculated  amplitude  and 
location  of  the  head  -  and  tail  waves  of  the  left-running  rarefaction,  the 
contact  surface  moving  to  the  right  and  the  shock  wave  moving  at  supersonic 
velocity  to  the  right  respectively.  The  solid  lines  are  the  solutions 
obtained  by  the  RCM  at  different  time  levels;  the  zero  numerical  diffusion 
feature  of  the  method  is  evident  in  the  'infinite'  resolution  of  the  contact 
discontinuity  and  the  shock,  and  the  dispersion  (phase  error)  is  within  one 
grid  spacing.  The  constant  states  are  perfectly  realized. 

It  is  these  features  of  the  method  that  make  it  very  attractive  for 
application  to  wave  rotor  type  flows,  since  the  successful  design  of  the 
device  is  predicated  on  being  able  to  accurately  compute  wave  arrival  times  at 
the  various  ports. 

2.3.2.  Wave  Turbine  Experiment 

Ref.  (10)  describes  the  wave  rotor  experimental  set  up  at  the 
Turbopropulsion  Laboratory.  Initial  tests  being  carried  out  currently  are 
with  the  wave  rotor  in  a  turbine  mode,  i.e.,  one  side  of  the  rotor  is  blocked 
off,  and  high  pressure  air  is  brought  onto  the  rotor  and  taken  off  again  from 


the  other  side.  The  passages  of  the  rotor  being  angled  at  60“  to  the  axis, 
the  180“  reversal  in  the  direction  of  the  fluid  flow  creates  an  angular 
momentum  change,  in  turn  generating  large  turbomachinery  work  coefficients. 
Fig.  (6)  shows  the  wave  diagram  computed  using  the  code.  The  movement  of  the 
rotor  is  from  top  to  bottom.  At  t=0  ,  the  high  pressure  air  is  brought 

into  contact  with  quiescent  atmospheric  air  in  the  rotor  passages,  at  point  a. 
This  corresponds  to  the  'piston'  inflow  boundary  condition  described  in 
section  2.2.3..  A  shock,  S  ,  is  generated  immediately,  (idealized  case  of 
instantaneous  cell  opening),  which  travels  from  the  right  to  the  left,  and 
strikes  the  solid  wall  at  the  left  end.  The  reflection  of  the  shock  takes 
place  at  point  b  according  to  the  solid  wall  boundary  condition  described  in 
section  2.2.1..  Behind  this  shock,  and  moving  at  a  slower  velocity  is  the 
contact  surface,  I  ,  which  penetrates  into  the  passage  only  a  fractional 
distance  before  encountering  the  reflected  shock,  RS,  at  point  c.  The 
reflected  shock  is  transmitted  through  the  contact  surface,  (bringing  the  flow 
to  a  near  halt),  and  reaches  the  right  side  at  point  d,  whereupon  the  inlet 
port  is  closed.  The  air  trapped  in  the  rotor  passages  is  now  at  a  high 
pressure  and  in  a  quiescent  state.  When  this  air  is  released  at  point  e  to  a 
low  pressure  region,  a  rarefaction  wave  is  generated,  R  ,  which  travels  to  the 
left,  spreading  out  in  the  process.  It  interacts  with  the  stationary  contact 
surface,  I  ,  setting  it  into  motion  again,  and  reflects  off  the  solid  wall  at 
the  left  as  RR  .  The  boundary  condition  imposed  at  point  e  is  the  outflow 
at  constant  static  pressure  condition  described  in  section  2.2.2..  The  outlet 
port  is  closed  at  a  time  when  the  exit  velocity  falls  to  about  half  its 
initial  value. 

This  experiment  embodies  two  fundamental  processes  in  wave  rotors: 
those  of  cell  filling  and  cell  emptying.  Almost  all  the  other  processes 
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typical  to  wave  rotors  are  combinations  of  the  cell  filling  and  cell  emptying 
unit  processes.  Comparison  of  the  ideal  computed  numbers  obtained  here  with 
experimental  data  will  provide  information  helpful  in  the  identification  and 
sources  of  losses. 

The  program  is  set  up  to  start  at  t*0  in  this  case,  with  initial 
data  provided  along  the  entire  passage,  i.e.,  from  x*0  to  x^O. 1863m  (the 
actual  length  of  the  wave  rotor  being  tested).  Since  the  passages  have 
quiescent  atmospheric  air  in  them  at  t>0  ,  the  initial  data,  of  course, 

describes  these  conditions.  Switches  for  the  left  and  right  boundaries 
describe  what  type  of  boundary  conditions  prevail  and  direct  the  program  to 
the  appropriate  subroutines.  These  switches,  designated  SWL  and  SWR  ,  for 
left  and  right  respectively,  are  assigned  integer  number  values  which 
correspond  to  the  numeric  value  of  the  particular  boundary  condition  they 
represent.  Thus,  if  the  left  boundary  is  a  solid  wall,  SWL*1  , 
corresponding  to  the  boundary  condition  subroutine  BCLl  .  In  this  example 
then,  the  initial  switch  settings  at  t-0  are  SWL*1  and  SWR«3  , 
corresponding  to  a  solid  wall  at  the  left  and  a  ’piston’  Inflow  at  the  right 
(which  starts  at  t-0  at  point  a).  At  point  d  ,  the  switches  are  reset  to 
SWL*1  and  SWR*1  due  to  the  closure  of  the  inlet  port.  At  point  e  ,  the 
switches  are  SWL”1  ,  SWR=2  ,  signifying  opening  of  the  exhaust  port  with 
outflow  at  a  constant  static  pressure.  The  whole  wave  diagram  can  thus  be 
packaged  into  a  ’module’  subroutine  and  called  from  the  main  program  with  a 
single  call  statement.  This  type  of  modularity  allows  for  wave  diagrams  of 
different  ’families’  to  be  developed  by  simply  calling  the  right  ’module’ 
subroutine. 

The  next  two  examples  illustrate  this  concept  as  they  deal  with  two 
very  different  types  of  wave  diagrams. 
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2.3.3.  General  Electric  Wave  Engine 


Fig.  (7)  shows  a  schematic  of  the  wave  diagram  constructed  for  the 
G.E.  wave  engine.  Briefly,  the  device  is  configured  for  a  gas  generator  mode 
of  operation,  with  counter flow  scavenging,  and  is  capable  of  producing  net 
shaft  power.  For  a  fuller  description  of  the  machine,  see  Ref.  (11).  In  this 
example,  fresh  charge  (air)  is  induced  into  the  rotor  (from  an  external 
reservoir)  through  the  wave  action  of  the  rarefaction  fan  originating  at  the 
exhaust  port  opening.  The  usefulness  of  the  rotor  is  gauged  by  the  net 
pressure  rise  across  the  machine,  i.e.,  the  ratio  of  the  total  exhaust 
pressure  to  total  (fresh  air)  inlet  pressure. 

For  performance  estimation  purposes,  it  is  sufficient  to  investigate 
only  the  exhaust  and  induction  processes  as  shown  in  Fig.  (8).  The  initial 
data  specified  is  as  follows:  the  exhausting  pressure  ratio  Pe/Po  »  the 
total  pressure  ratio  across  the  rotor  Pte^Pta  assumed  total 

temperature  ratio  T^g/T^a  .  In  this  particular  cycle,  the  amount  of  fresh 
charge  induced  in  is  ideally  equal  to  the  gases  exhausted  out,  i.e., 

®in  *  “out  »  this  mass  balance  is  carried  out  after  each  computation  to 

correct  the  assumed  temperature  ratio  T^g/Tta  (which  otherwise  constitutes 
overspecification  of  the  initial  conditions). 

The  calculation  starts  at  t=0  ,  with  initial  data  consistent  with 
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the  chosen  pressure  and  temperature  ratios  specified  along  the  passage  length. 
Initial  switch  settings  are  SWL=1  and  SWR=2  for  the  solid  wall  boundary  at 
the  left  and  the  exhaust  to  a  constant  pressure  at  the  right.  As  shown  in  the 
figure,  a  rarefaction  fan  is  generated,  propagating  to  the  left  and  reflecting 
off  the  solid  wall.  At  time  t=T^  ^  the  pressure  at  the  wall  has  been 
reduced  to  that  outside  the  passage,  pta  »  which  is  when  the  inlet  port  is 
opened.  The  switches  are  now  set  to  SWL*4  and  SWR=2  for  isentropic  inflow 


from  an  external  reservoir  at  the  left,  and  still  outflow  at  a  constant 
pressure  at  the  right.  The  exhaust  port  is  closed  at  time  t=T2  which 
corresponds  to  the  exit  velocity  having  dropped  off  to  approximately  half  its 
steady  state  value  at  the  beginning  of  the  exhaust  process.  Now  the  switches 
are  set  to  SWL*4  and  SWR=1  ,  for  the  solid  wall  condition  at  the  right. 

The  sudden  closure  of  the  exhaust  port  generates  a  'hammer'  shock  travelling 
to  the  left,  interacting  with  the  incoming  interface  (shown  by  dashed  line), 
and  reaching  the  passage  end  at  t=T4  at  which  time  the  inlet  port  is  closed, 
with  the  switches  being  reset  to  SWL=1  and  SWR=1  .  Note  the  reflected 
shock  travelling  from  left  to  right  generated  at  the  interaction  of  the 
contact  surface  and  the  hammer  shock. 

Once  this  solution  is  obtained,  integration  of  the  mass  flux  through 
the  inlet  and  exhaust  ports  is  carried  out  and  if  the  two  numbers  do  not 
match,  the  assumed  temperature  ratio  T^e/Tta  is  adjusted  in  the  initial 
data,  till  such  time  as  m^n  =  %)ut  ♦ 

This  calculation  is  sufficient  for  performance  analyses:  if  the 
entire  wave  diagram  has  to  be  worked  out,  then  at  a  time  t  >  T3  ,  hot  gas 
from  the  combustion  chamber  is  brought  onto  the  rotor  (the  boundary  condition 
corresponding  to  'piston'  inflow)  on  the  right  hand  side.  This  would 
generate  the  shock  to  compress  the  induced  air  and  when  this  shock  reached 
Che  left  end,  the  transfer  port  (see  Fig.  7)  would  be  opened  for  such  time  it 
takes  for  the  compressed  air  to  be  completely  scavenged  out  of  the  rotor. 

Fig.  (9)  shows  some  performance  curves  obtained  using  the  procedure  outlined 
above.  In  Figs.  (10a,  b,  c)  are  shown  three  sets  of  flow  parameters  at 
different  time  steps  corresponding  to  the  inlet  port  just  opening,  the 
exhaust  port  closing  and  the  inlet  port  closing;  the  qualitative  distributions 
of  the  flow  parameters  in  the  passage  are  immediately  seen  to  be  accurate  when 


compared  with  the  wave  diagram  shown  in  Fig.  (8).  Of  interest  is  the  set  of 
plots  for  the  time  step  when  the  inlet  port  has  just  been  closed.  The  flow 
between  the  end  of  the  passage  and  the  location  of  the  interface  is  seen  to 
be  quite  non-uniform  in  the  density  plot.  At  the  same  time,  the  shock 
reflected  from  the  interface  has  reached  the  right  side  and  reflected  off  the 
solid  wall.  These  considerations  help  to  decide  optimum  port  opening  and 
closing  times.  For  example.  Fig.  (11)  shows  what  happens  if  the  inlet  port  is 
not  closed  at  just  the  time  the  shock  reaches  the  end,  but  rather  at  some 
short  time  later.  The  shock  now  sees  an  open  boundary  and  reflects  off  as  an 
expansion  to  match  the  high  pressure  behind  it  with  the  incoming  total 
pressure  which  is  at  a  lower  value.  This  reflected  expansion  is  manifested  in 
the  pressure,  density  and  velocity  plots  of  the  figure. 

The  entire  sequence  of  wave  interactions  of  this  example  is  computed 
by  the  RCM  without  the  implementation  of  artificial  viscosity  or  artificial 
compression  methods,  or  tracking  and  capturing  schemes.  This  'hands  off' 
feature  of  the  method  renders  it  eminently  useful  for  fast  preliminary 
evaluations  of  complex  wave  diagrams  for  the  application  at  hand. 

The  next  example  computes  an  idealized  wave  diagram  for  the  nine-port 
pressure  exchanger  concept  proposed  by  Spectra  Technology,  Ref.  (12). 

2.3.4.  Spectra  Technology  Pressure  Exchanger 

Fig.  (4)  shows  the  ideal  wave  diagram  for  the  nine-port  pressure 
exchanger.  This  configuration  is  a  good  case  example  to  compute  with  the  RCM 
because  of  the  different  types  of  boundary  conditions  that  need  to  be  dealt 
with  in  the  evaluation  of  the  cycle.  Th<^  computation  is  started  at  t=0  ,  at 

the  point  of  high  pressure  hot  gas  inlet  (driver  gas  inlet).  In  the  manner 
described  in  the  G.E.  wave  engine  example,  the  initial  data  is  prescribed  for 


U 


I  . 


Li 
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the  entire  passage  at  this  time  step  and  the  boundary  condition  switches  are 
initially  set  at  SWL=1  and  SWR^S  for  the  solid  wall  at  the  left,  and  the 
'piston'  inflow  at  the  right  hand  end.  Since  there  is  a  multiplicity  of  types 
of  boundary  conditions,  e.g. ,  three  outflow  ports,  an  index,  JCOUNT,  is  used 
to  ensure  proper  sequencing  of  the  switches.  The  following  table  is 
presented  as  an  example  of  the  settings  of  the  switches  to  carry  out 
calculations  for  one  cycle.  The  inflow  and  outflow  port  conditions  are  those 


proposed  by  Spectra  Technology  for  their  idealized  diagram. 


TIME  STEP,  N 

JCOUNT 

SML 

SWR 

REMARKS 

0 

0 

1 

3 

CYCLE  STARTS.  HP  GAS  INLET  PORT  OPENS 

500 

1 

2 

3 

HP  AIR  OUTLET  PORT  OPENS 

1408 

2 

2 

1 

HP  GAS  INLET  PORT  CLOSES 

1765 

3 

5 

1 

HP  AIR  OUTLET  PORT  CLOSES.  TUNING  PORT 
LI  OPENS 

1816 

4 

2 

1 

TUNING  PORT  LI  CLOSES.  IP  GAS  OUTLET 
PORT  OPENS  (PORT  El) 

2069 

5 

2 

5 

TUNING  PORT  R1  OPENS 

2261 

6 

2 

1 

TUNING  PORT  R1  aOSES 

2595 

7 

5 

1 

IP  GAS  OUTLET  PORT  CLOSES.  TUNING  PORT 
L2  OPENS 

2636 

8 

2 

1 

TUNING  PORT  L2  CLOSES.  LP  GAS  OUTLET 
PORT  OPENS  (PORT  E2) 

3029 

9 

2 

5 

TUNING  PORT  R2  OPENS 

3237 

10 

2 

m 

TUNING  PORT  R2  CLOSES.  LP  AIR  INLET 
PORT  OPENS 

4961 

11 

1 

H 

LP  GAS  OUTLET  PORT  CLOSES 

5529,0 

0 

1 

3 

LP  AIR  INLET  PORT  CLOSES.  CYCLE 
COMPLETED 

The  total  cycle  time  as  calculated  by  the  RCM  is  3.0676  mseconds,  which 


compares  well  with  the  time  computed  by  Spectra  Technology  (using  the 
FCT-SHASTA  algorithm)  of  3.07  mseconds.  The  execution  time  on  an  IBM 
370-3033AP  for  the  5529  steps  computed  in  the  example  above  was  3  minutes 
38  seconds,  including  the  I/O  operations  and  the  graphics. 


Figs.  (12a,  b,  c)  show  three  sets  of  plots  of  the  flow  parameters  for 
the  following  cases:  a)  the  H.P.  air  outlet  port  opens  on  time,  l.e.,  just  as 

the  shock  reaches  the  left  end  of  the  passage,  b)  the  port  opens  before  the 

shock  has  reached  the  end,  and  c)  the  port  opens  after  the  shock  has  reached 

the  end.  The  constant  pressure  and  velocity  states  that  prevail  in  the 
passage  just  after  the  shock  has  reached  the  left  end  (time  'section'  line 
in  wave  diagram),  are  perfectly  realized  in  Fig.  (12a),  while  the  contact 
surface  is  at  the  location  shown  by  the  sharp  discontinuities  in  the  density 
and  entropy  plots.  Should  the  inlet  port  be  opened  earlier,  e.g.,  at  the  time 
level  shown  by  in  the  wave  diagram,  what  happens  is  as  follows:  the 

pressure  in  the  passage  is  still  at  the  pre~compressed  level  and  this  comes 
into  contact  with  the  pressure  level  in  the  port  which  is  considerably  higher, 
resulting  in  a  shock  propagating  into  the  passage,  colliding  with  the  left 
moving  shock  and  raising  the  overall  pressure  level  to  -'S.O  as  shown  in 
Fig.  (I2b).  However,  as  soon  as  the  left  moving  shock  reaches  the  end,  it  now 
encounters  an  open  boundary  with  conditions  that  do  not  match  those  behind  the 
shock,  resulting  in  a  rarefaction  fan  being  generated,  which  propagates  to  the 
right.  This  expansion  fan,  travelling  at  sonic  velocity  relative  to  the  gas 
into  which  it  is  propagating,  soon  overtakes  the  right  moving  shock  which  is 
travelling  at  a  subsonic  velocity  relative  to  the  same  gas.  This  interaction 
results  in  an  attenuation  of  both  the  rarefaction  as  well  as  the  shock  wave. 
Note  that  the  overall  pressure  and  velocity  levels  behind  the  rarefaction  are 
about  the  same  as  for  case  a),  i.e.,  the  effects  of  the  mismatch  are  not  very 
significant  at  the  outlet  port.  However,  should  the  right  moving  pressure 
perturbations  of  case  b)  not  attenuate  each  other  significantly  before  they 
reach  the  right  hand  end,  the  consequences  could  be  severe  for  the  overall 
wave  diagram,  since  this  will  lead  to  further  (unwanted)  wave  reflections. 
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Fig,  (12c)  shows  what  occurs  if  the  outlet  port  is  opened  too  late, 
corresponding  to  time  level  1^+  on  the  wave  diagram.  Now  the  left 
travelling  shock  encounters  a  wall  boundary  condition  on  reaching  the  left 
end  and  reflects  off  as  a  shock,  effectively  doubling  the  pressure  level 
behind  it  (  >3.5  in  pressure  plot  of  Fig,  (12c)).  When  the  outlet  port  opens, 
there  is  again  a  mismatch  of  conditions  in  the  port  and  in  the  passage,  with 
the  pressure  level  in  the  passage  being  considerably  higher  than  that 
prescribed  for  the  outlet  port.  A  rarefaction  wave  is  generated  which 
propagates  to  the  right  and  overtakes  the  reflected  shock.  The  same  criterion 
holds  for  this  case  too,  l.e.,  the  ensuing  attenuation  of  these  pressure 
pulses  should  occur  before  they  reach  the  right  hand  end,  preferably  even 
before  they  reach  the  interface  still  propagating  towards  the  left  at  the 
flow  velocity. 

The  considerations  above  give  a  preview  of  the  nature  of  decisions 
required  in  the  successful  design  of  a  wave  rotor  device.  It  is  clear  that 
quite  a  few  iterations  are  involved  in  the  process  of  designing  a  viable  wave 
diagram  for  a  particular  application,  and  each  iteration  entails  calculating 
two  or  more  complete  cycles  to  ensure  ’closure’  or  repeatability  of  the  cycle. 
A  fast  solver  like  the  RCM  allows  reaching  an  idealized  ’base’  design  quickly 
and  inexpensively. 

Appendix  A  is  a  listing  of  the  program  in  its  present  development 
stage.  As  mentioned  earlier,  no  attempt  has  been  made  to  optimize  the 
program,  either  for  storage  requirements  or  for  execution. 

Appendix  B  gives  a  description  of  the  structure  of  the  program,  a 
listing  of  the  important  variables,  the  subroutines  and  the  function 
subprograms.  A  step  by  step  guide  is  also  included  to  set  up  and  run  the 


program 


3.  DISCUSSION  AND  RECOMMENDATIONS 


3.1.  Discussion 

For  meeting  the  criteria  listed  In  the  Introduction,  In  one  dimension, 
Glimm's  method  or  the  RCM  appears  to  be  superior  to  any  difference  method. 

For  wave  rotor  type  applications,  where  discontinuities  need  to  be  computed 
with  sharpness,  the  'Infinite'  resolution  of  such  discontinuities  Inherent  In 
the  RCM  make  It  a  natural  choice  to  carry  out  Ideal  flow  calculations  for 
preliminary  design  purposes.  Boundary  conditions  can  be  Implemented  quite 
easily  and  do  not  require  Information  from  points  outside  the  domain  of 
dependance  as  Is  the  case  In  some  finite  difference  schemes.  The  van  der 
Corput  sampling  technique  results  In  the  best  possible  representation  of  the 
wave  propagation,  which  Is  essential  for  the  correct  representation  of 
continuous  waves,  particularly  those  produced  by  nonlinear  Interactions. 

The  method,  however,  is  not  recommended  to  solve  for  flows  with  real 
effects  such  as  friction,  heat  transfer  and  area  change,  or  to  be  extended  to 
multl-dlmenslonal  flows.  Although  considerable  research  is  being  done  to 
rigorously  extend  the  method  to  such  flows,  with  some  degree  of  success  (see 
Refs.  1,  4,  13),  the  present  state  of  development  is  not  mature  enough  to 
ensure  a  useful  practical  code  as  the  outcome. 

3.2.  Recommendations 

Many  options  are  available  for  one  wishing  to  develop  either  a  1-D  code 
with  real  effects  and/or  a  multl-dlmenslonal  code  for  wave  rotor  type 
applications.  The  author  prefers  to  recommend  numerical  formulations  which 
are  dependent  on  the  solution  of  Rlemann  problems,  such  as  the  Godunov  method; 
the  motivating  reason  for  this  preference  Is  that  a  Rlemann  problem 
constitutes  the  solution  of  a  discontinuity  in  the  flow  in  terms  of  other 


discontinuities  (if  any  are,  indeed,  present),  and  the  scheme  is  thus 
intrinsically  suited  for  solving  such  flows;  on  the  other  hand,  the  other 
schemes,  in  general,  require  to  be  made  aware  of  discontinuities  in  the  flow 
through  some  external  device,  and  then  treat  them  through  other  artificial 
devices. 

A  second-order,  quasi  one-dimensional  (variable  cross-sectional  area) 
scheme  has  recently  been  developed  by  Ben-Artzi  and  Falcovitz  (Ref.  14).  The 
method  is  based  on  the  exact  solution  of  'generalized  Rlemann  propblems',  and 
has  demonstrated  very  good  results;  it's  least  accurate  approximation  is 
equivalent  to  Godunov's  first  order  method  (Ref.  9).  The  resolution  of  shocks 
and  other  dlsconltinultles  and  singularities  of  the  flow  field  is  also  high. 
Extension  to  more  than  one  dimension  appears  to  be  straightforward  through  the 
use  of  operator  splitting  techniques,  but  has  as  yet  not  been  tried 


extensively 
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Fig. 12  a  :  Distribution  of  Flow  Parameters  in 
Rotor  Passage  when  the  H.P.  Air 
Outlet  Port  Opens  on  Time. 
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PROGRAM  RCM  WITH  VAN  DER  CORPUT  SAMPLING  AND  SINGLE  TIME  STEP 

RCM00030 

'■ 

INTEGER  QPRINT,  QSTOP,  SWL,  SWR 

RCM00040 

DIMENSION  XX(6),YY(6) 

RCM00050 

i'.A 

DIMENSION  XARRAY(IOO) 

DIMENSION  WN0RM(12) ,IDIGT(12) 

RCM00060 

RCM00070 

V-; 

DIMENSION  P(203),R(203),U(203) ,A(203) ,S(203),X(203) 

RCM00080 

S 

COMMON/ SUBS/P,R,U, A, S,X 

COMMON / GLIMMl / PGLIM , RGLIM , UGLIM , PL , RL , UL , PR , RR , UR , AL , AR , GL , GR , EPS 

RCM00090 

RCMOOlOO 

COMMON/ GLIMM2 / DT , DX , XI 

RCMOOllO 

p? 

COMMON/ FUNl/ G , PA , RA , UA , RB , RMU 

RCM00120 

COMMON / SAMPLE / WNORM , IDIGT 

RCM00130 

COMMON  XARRAY.Nl 

RCM00140 

CALL  COMPRS 

RCM00150 

CALL  BLOWUP (0.5) 

CALL  PAGE (11. 0,8. 5) 

CALL  HWSCAL( ' SCREEN' ) 

RCM00160 

RCM00170 

RCM00180 

DATA  K,SWL,SWR/500,1,3/ 

RCM00190 

DATA  N,CFLNUM,TTOTAL/0, 0.60,0.0/ 

RCM00200 

DATA  PSEXIT,PSINL,PSINR,RINL,RINR/ 116954. ,3770000. ,3819952.50,6.8 

,RCM00210 

f. 

*6.800/ 

RCM00220 

DATA  PSOUTl , PSOUT2 , PSOUT3/ 38 19952 . 5,2431800 . 0 , 1530007 . 5/ 

DATA  PTOTIN , RTOTIN/ 1656663.8,7.486/ 

RCM00230 

RCM00240 

ts 

DATA  PREF,RREF,XREF/ 1656663. 8, 7. 486, 0.1800/ 

RCM00250 

I* ' 

G=1.4 

GL=1.4 

RCM00260 

RCM00270 

GRsl.4 

RCM00280 

‘ 

EPS=l.E-06 

RCM00290 

N 

QSTOP=20 

NlsO 

RCM00300 

RCM00310 

JCOUNT=0 

RCM00320 

• 

KC0UNT=0 

UEXMAX=0. 

DX=0.01 

RCM00330 

RCM00340 

RCM00350 

AREF=SQRT(PREF/RREF) 

RCM00360 

TIMREF=XREF/AREF 

RCM00370 

RMU=SQRT((G-1. )/(G+l. )) 

RCM00380 

X(1)=-0.5*DX 

RCM00390 

ZETA=WDP(1) 

XIl=DX*(WDP(0)-0.5) 

RCM00400 

RCM00410 

DO  25  1=2,203 

RCM00420 

X(I)=X(I-1)+0.5*DX 

25  CONTINUE 

RCM00430 

RCM00440 

- 

DO  35  1=1,100 

RCM00450 

.  ^ 

XARRAY(I)=X(I*2+1) 

RCM00460 

35  CONTINUE 

INITIAL  DATA 

RCM00470 

RCM00480 

CALL  INITl 

RCM00490 

CALL  INIT2L(PSEXIT) 

CALL  INIT2R(PSEXIT,PREF,RREF) 

RCM00500 

RCM00510 

'*• 

CALL  INIT3L(PSINL,RINL) 

RCM00520 

CALL  INIT3R(PSINR,RINR) 

NONDIMENS lONALI ZATION 

RCM00530 

RCM00540 

DO  30  1=1,203,2 

RCM00550 

P(1)=P(I)/PREF 

RCM00560 

R(I)=R(I)/RREF 
U(I)=U(I)/AREF 
A(I)=A(I)/AREF 
S(I)=ALOG(P(I)/R(I)**G) 
30  CONTINUE 

CALL  PLOTl(K) 

DO  40  J=1,K 
N=N+1 

XII=DX*(WDP(0)-0.5) 

QPRINT=N/50 


RCM00570 

RCM00580 

RCM00590 

RCM00600 

RCM00610 

RCM00620 

RCM00630 

RCM00640 

RCM00650 

RCM00660 


*  •  "j 


DT=100. 


RCM00670 


DO  50  1=1,203,2 

DTT=CFLNUM*DX/(2.*AMAX1(ABS(A(I)+U(I)) ,ABS(A(I)-U(I)))) 
DT=AMIN1(DTT,DT) 

50  CONTINUE 

TTOTAL=TTOTAL+DT 

TIME=TTOTAL*TIMREF 

XI=-XII 


RCM00680  V- 
RCM00690  --I: 

RCM00700  i::!: 

RCM00710  P 
RCM00720  I-V 
RCM00730 
RCM00740  :•  V 


DO  60  1=1,201,2 
PL=P(I) 

RL=R(I) 

UL=U(I) 

PR=P(I+2) 

RR=R(I+2) 

UR=U(I+2) 

XITEMP=XI 

IF(I.EQ.l)  XI=ABS(XI) 

IF((I.EQ.201).AND. (XI.GT.0.0))  XI=-XI 
CALL  GLIMM ( QSTOP , P STAR , USTAR , ASTAR ) 

XI=XITEMP 
P(I*1)=PGLIM 
R(I+1)=RGLIM 
U(I+1)=UGLIM 
60  CONTINUE 

DO  70  1=1,201,2 
IF(XI.LT.O.)  GOTO  80 
P(I+2)=P(I+1) 

R(I>2)=R(I+1) 

U(I+2)=U(I+1) 

A(I+2)=SQRT(G*P(I+2)/R(I+2)) 

S(I+2)=AL0G(P(I+2)/R(I+2)**G) 

GOTO  70 
80  P(I)=P(I+1) 

R(I)=R(I>1) 

U(I)=U(I+1) 

A(I)=SQRT(G*P(I)/R(I)) 

S(I)=ALOG(P(I)/R(I)**G) 

70  CONTINUE 

CALL  GE ( SWL , SWR , N , TTOTAL , TIME , UEXMAX , PTOTIN , PREF ) 

CALL  DETON ( SWL , SWR , N , QPRINT , TTOTAL , TIME ) 

CALL  SPCTRA(N , SWL , SWR , TIME , UEXMAX , PSEXIT , PSOUTl , PSOUT2 , PSOUT3 , J 
*COUNT , QPRINT , TTOTAL ,KCOUNT ) 

IF(SWL.EQ.  1',  CALL  BCLl 
IF(SWL.EQ.2)  CALL  BCL2 (PSEXIT , PREF) 


RCM00750 

RCM00760 

RCM00770 

RCM00780 

RCM00790 

RCM00800 

RCM00810 

RCM00820 

RCM00830 

RCM00840 

RCM00850 

RCM00860 

RCM00870 

RCM00880 

RCM00890 

RCM00900 

RCM00910 

RCM00920 

RCM00930 

RCM00940 

RCM00950 

RCM00960 

RCM00970 

RCM00980 

RCM00990 

RCMOIOOO 

RCMOlOlO 

RCM01020 

RCM01030 

RCM01040 

RCM01050 

RCM01060 

RCM01070 

RCM01080 

RCM01090 

RCMOllOO 


IF(SWL.EQ.3)  CALL  BCL3(PSINL,RINL,PREF,RREF) 

IF(SWL.EQ.4)  CALL  BCL4 (PTOTIN .RTOTIN ,PREF ,RREF) 

IF(SWL.EQ.5)  CALL  BCL5 
IF(SWR.EQ.l)  CALL  BCRl 
IF(SWR.EQ.2)  CALL  BCR2 (PSEXIT , PREF) 

IF(SWR.EQ.3)  CALL  BCR3 (PSINR,RINR,PREF,RREF) 

IF(SWR.EQ.4)  CALL  BCR4( PTOTIN, RTOTIN, PREF, RREF) 

IF(SWR.EQ.5)  CALL  BCR5 

IF  ((N.EQ. (50*QPRINT)) .AND. (N.GE.O))  CALL  PLOT2(N,K) 

40  CONTINUE 

CALL  ENDPL(O) 

CALL  DONEPL 

STOP 

END 

SUBROUTINE  GLIMM ( QSTOP , PSTAR , USTAR , ASTAR ) 

INTEGER  Q, QSTOP 
REAL  ML,MR,MLN,MRN 

COMMON / GLIMMl / PGLIM , RGLIM , UGLIM , PL , RL , UL , PR , RR , UR , AL , AR , GL , GR , EPS 
COMMON / GLIMM2 / DT , DX , XI 
DATA  Q,ML,MR/0, 100. , 100. / 

PSTAR=0.5*(PL+PR) 

COEFL=  SQRT ( PL*RL ) 

COEFR=  SQRT ( PR*RR ) 

ALPHA=1, 

BEGIN  GODUNOV  ITERATION 
30  Q=Q+1 

IF(PSTAR.LT.EPS)  PSTAR=EPS 
COMPUTE  NEXT  ITERATION  FOR  ML  AND  MR 
MLN= COEFL*PHI (PSTAR , PL ) 

MRN=  C0EFR*PHI (PSTAR , PR ) 

DIFML=ABS(MLN-ML) 

DirMR=ABS(MRN-MR) 

ML=MLN 

MR=MRN 

COMPUTE  NEW  PSTAR 
PTIL=  PSTAR 

PSTAR= (UL-UR+PL/ML+PR/MR) / ( 1 . /ML+ 1 . /MR) 

PSTAR=ALPHA*PSTAR+ ( 1 . - ALPHA )*PTIL 
IF (Q.LE. QSTOP)  GOTO  10 
IfUbS(PSTAR-PTIL)  .LT.EPS)  GOTO  20 
COMPUTE  NEW  ALPHA 
ALPHA=0 . 5*ALPHA 
Q=0 

IF( (1. -ALPHA) .LT. EPS)  GOTO  20 
10  IF (DIFML.GE.EPS)  GOTO  30 
IF (DIFMR.GE.EPS)  GOTO  30 
END  OF  GODUNOV  ITERATION;  COMPUTE  USTAR 
20  USTAR= (PL-PR+ML*UL+MR*UR) / (ML+MR) 

BEGIN  SAMPLING  PROCEDURE 
IF  (XI.LT.USTAR*DT)  GO  TO  40 

RIGHT  SIDE;  SELECT  CASE  OF  SHOCK  OR  EXPANSION 
IF  (PSTAR.lt. PR)  GO  TO  50 
RIGHT  WAVE  IS  A  SHOCK  WAVE 
WR=UR+MR/RR 
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IF  (XI.LT.WR*DT)  GO  TO  60 

RIGHT  OF  RIGHT  SHOCK  CASE 

RGLIM=RR 

PGLIM=PR 

UGLIM=UR 

RETURN 

LEFT  OF  RIGHT  SHOCK  CASE 
60  RGLIM=-MR/(USTAR-WR) 

PGLIM=PSTAR 

UGLIM=USTAR 

RETURN 

RIGHT  WAVE  IS  A  RAREFACTION  WAVE 
50  CONST=PR/RR**GR 

RSTAR= ( PSTAR/ CONST ) ** ( 1 . / GR ) 

ASTAR=  SQRT (GR*PSTAR/RSTAR) 

AR=SQRT(GR*PR/RR) 

IF  (XI.GE. (USTAR+ASTAR)*DT)  GO  TO  70 

LEFT  OF  RIGHT  FAN  CASE 

RGLIM=RSTAR 

UGLIM=USTAR 

PGLIM= PSTAR 

RETURN 

SELECT  RIGHT  OF  FAN  OR  IN  FAN 
70  IF  (XI.GE. (UR+AR)*DT)  GO  TO  80 
IN  RIGHT  FAN  CASE 

UGLIM=2./(GR+1. )*(XI/DT-AR+(GR-1. )/2.*UR) 

RGLIM= ( ( ( AR+ (GR- 1 . ) / 2 . * (UGLIM-UR) ) **2 . ) / (GR*CONST ) ) ** ( 1 . / (GR- 1 , 

PGLIM=CONST*RGLIM**GR 

RETURN 

RIGHT  OF  RIGHT  FAN  CASE 
80  RGLIM=RR 
PGLIM=PR 
UGLIM=UR 
RETURN 

LEFT  SIDE;  SELECT  CASE  OF  SHOCK  OR  RAREFACTION 
40  IF  (PSTAR.LT.pl)  GO  TO  90 
LEFT  WAVE  IS  A  SHOCK  WAVE 
WL=UL-ML/RL 

IF  (XI.GE.WL*DT)  GO  TO  100 

LEFT  OF  LEFT  SHOCK  CASE 

RGLIM=RL 

PGLIM=PL 

UGLIM=UL 

RETURN 

RIGHT  OF  LEFT  SHOCK  CASE 
100  RGLIM=ML/(USTAR-WL) 

PGLIM= PSTAR 
UGLIM=USTAR 
RETURN 

LEFT  WAVE  IS  A  RAREFACTION  WAVE 
90  CONST = PL/ RL**GL 

RSTAR= (PSTAR/ CONST )**( 1. /GL) 

ASTAR= SQRT (GL*P STAR /RSTAR) 

AL=SQRT(GL*PL/RL) 
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IF  (XI.LT. (USTAR-ASTAR)*DT)  GO  TO  110 

RIGHT  OF  LEFT  FAN  CASE 

RGLIM=RSTAR 

PGLIM=PSTAR 

UGLIM=USTAR 

RETURN 

SELECT  LEFT  OF  FAN  OR  IN  FAN  CASE 
110  IF  (XI.LT. (UL-AL)*DT)  GO  TO  120 
IN  LEFT  FAN  CASE 

UGLIM=2./(GL+1. )*(AL+(GL-1. )/2.*UL+XI/DT) 

RGLIM= ( ( (AL+ (GL-1. )/2 . *(UL-UGLIM) )**2 . ) / (GL*C0NST) )**( 1 . / (GL-1. 

PGLIM=CONST*RGLIM**GL 

RETURN 

LEFT  OF  LEFT  FAN  CASE 
120  RGLIM=RL 
PGLIM=PL 
UGLIM=UL 
RETURN 
END 

FUNCTION  PHI(Y,Z) 

REAL  RMU 

COMMON/ FUNl/ G , PA , RA , UA , RB , RMU 

EPS=l.E-06 

PARAM=Y/Z 

IF  (ABS(l.-PARAM).GE.EPS)  GO  TO  10 
PHI=SQRT(G) 

RETURN 

10  IF  (PARAM.GE.l. )  GO  TO  20 

PHI=(G-1. )/2.*(l.-PARAM)/(SQRT(G)*(l.-PARAM**((G-l. )/(2.*G)))) 
RETURN 

20  PHI=SQRT((G+1. ) /2 . *PARAM+ (G- 1 . )/2. ) 

RETURN 

END 

FUNCTION  PHIl(PB) 

REAL  RMU 

COMMON/ FUNl/ G , PA , RA , UA , RB , RMU 

PHI 1 = ( PB - PA ) *SQRT ( ( 1 . - RMU**2 . ) / ( RA* ( PB  +  RMU**2 . *PA ) ) ) 

RETURN 

END 

FUNCTION  PSI(PB) 

REAL  RMU 

COMMON/ FUNl/ G , PA , RA ,UA , RB , RMU 

P  S I =  S  QRT ( 1 . - RMU**4 . ) / RMU**2 . / S  QRT ( RA ) *P A** (l./(2.*G))*(PB**((G- 
*/(2.*G))-PA**((G-l. )/(2.*G))) 

RETURN 

END 

SUBROUTINE  INITl 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 

COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 

COMMON/ SUBS/P, R,U, A, S,X 

DO  10  1=1, 9, 2 

P(I)=810600.00 

R(I)=0.7132 

U(I)=6A4.4 
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A(I)=SQRT(G*P(I)/R(I)) 

10  CONTINUE 

DO  20  1=11,203,2 

P(I)=101325.0 

R(I)=1.22 

uu)=o.o 

A(I)=SQRT(G*P(I)/R(I)) 

20  CONTINUE 
RETURN 
END 

SUBROUTINE  INIT2R ( PSEXIT , PREF , RREF ) 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 

COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 

COMMON/ SUBS/ P,R,U, A, S,X 

DO  10  1=3,201,2 

P(I)=PREF 

R(I)=RREF 

UU)  =  0-0 

A(I)=SQRT(P(I)*G/R(I)) 

10  CONTINUE 
P(1)=P(3) 

R(1)=R(3) 

U(l)=-U(3) 

A(1)=SQRT(G*P(1)/R(1)) 

P(203)=PSEXIT 

R(203)=R(201) 

PA=P(201) 

RA=R(201) 

UA=U(201) 

PB=P(203) 

RB=R(203) 

IF(PA.GT.PB)  GO  TO  20 
U(203)=UA-PHI1(PB) 

GO  TO  30 

20  U(203)=UA-PSI(PB) 

30  A(203)=SQRT(G*P(203)/R(203)) 

RETURN 

END 

SUBROUTINE  INIT2L(PSEXIT) 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 

COMMON / FUNl / G , PA , RA , UA , RB , RMU 

COMMON/ SUBS/P, R,U, A, S,X 

DO  10  1=3,201,2 

P(I)=285080.0 

R(I)=0.897 

U(I)=0.0 

A(I)=SQRT(G*P(I)/R(I)) 

10  CONTINUE 
RETURN 
END 

SUBROUTINE  INIT3L(PSINL,RINL) 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203),S(203) ,X(203) 
COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 
COMMON/ SUBS /P,R,U, A, S,X 
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DO  10  1=3,201,2 
P(I)=2390000.0 
R(I)=9.787 
U{I)=0.0 

A(I)=SQRT(G*P(I)/R(I)) 

10  CONTINUE 
P(1)=PSINL 
R(1)=RINL 
PA=P(3) 

RA=R(3) 

UA=U(3) 

PB=P(1) 

U(1)=UA+PHI1(PB) 

A(1)=SQRT(G*P(1)/R(1)) 

P(203)=P(201) 

R(203)=R(201) 

U(203)=-U(201) 

A(203)=SQRT(G*P(203)/R(203)) 

RETURN 

END 

SUBROUTINE  INIT3R(PSINR, RINR) 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 

COMMON/ SUBS/P,R,U, A, S,X 

COMMON/ FUNl/G , PA , RA , UA , RB ,  RMU 

DO  10  1=3,201,2 

P(I)=2421667.5 

R(I)=9.787 

U(I)=0.0 

A(I)=SQRT(G*P(I)/R(I)) 

10  CONTINUE 

P(203)=PSINR 

PA=P(201) 

RA=R(201) 

UA=U(201) 

PB=P(203) 

U(203)=UA-PHI1(PB) 

R(203)=RINR 

A(203)=SQRT(G*P(203)/R(203)) 

P(1)=P(3) 

R(1)=R(3) 

U(l)=-U(3) 

A(1)=SQRT(G*P(1)/R(1)) 

RETURN 

END 

SUBROUTINE  BCLl 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 
COMMON/ SUBS/P,R,U, A, S,X 
COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 
P(1)=P(3) 

R(1)=R(3) 

U(l)=-U(3) 

A(1)=SQRT(G*P(1)/R(1)) 

RETURN 

END 
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SUBROUTINE  BCRl 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 
COMMON/ SUBS/P,R,U, A, S,X 
COMMON/ FUNl/G , PA , RA ,UA , RB , RMU 
P(203)=P(201) 

R(203)=R(201) 

U(203)=-U(201) 

A(203)=SQRT(G*P(203)/R(203)) 

RETURN 

END 

SUBROUTINE  BCL2 (PSEXIT , PREF) 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 

COMMON/ SUBS/P, R,U, A, S,X 

COMMON/ FUNl/ G , PA , RA , UA , RB , RMU 

P(1)=PSEXIT/PREF 

R(1)=R(3) 

U(1)=U(3) 

A(1)=SQRT(G*P(1)/R(1)) 

RETURN 

END 

SUBROUTINE  BCR2 (PSEXIT , PREF) 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203) ,S(203),X(203) 

COMMON/ SUBS/P,R,U, A, S,X 

COMMON/ FUNl/G , PA , RA ,UA , RB , RMU 

P(203)=PSEXIT/PREF 

R(203)=R(201) 

PA=P(201) 

RA=RU01) 

UA=U(201) 

PB=P(203) 

RB=R(203) 

IF(PA.GT.PB)  GO  TO  10 
U(203)=UA-PHI1(PB) 

GO  TO  20 

10  U(203)=UA-PSI(PB) 

20  A(203)=SQRT(G*P(203)/R(203)) 

RETURN 

END 

SUBROUTINE  BCL3 ( PSINL , RINL , PREF , RREF ) 

DIMENSION  P(203)  ,R(203)  ,U(203)  ,A(203hS(203)  ,X(203) 

COMMON/ SUBS/P,R,U, A, S,X 

COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 

P(1)=PSINL/PREF 

R(1)=RINL/RREF 

PA=P(3) 

RA=R(3) 

UA=U(3) 

PB=P(1) 

U(1)=UA>PHI1(PB) 

A(1)=SQRT(G*P(1)/R(1)) 

RETURN 

END 

SUBROUTINE  BCR3 ( PSINR , RINR , PREF , RREF ) 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203),S(203),X(203) 
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COMMON/ SUBS/P,R,U, A, S,X 

COMMON/ FUNl / G , PA , RA , UA , RB ,  RMU 

P(203)=PSINR/PREF 

R(203)=RINR/RREF 

PA=P(201) 

RA=R(201) 

UA=U(201) 

PB=P(203) 

U(203)=UA-PHI1(PB) 

A(203)=SQRT(G*P(203)/R(203)) 

RETURN 

END 

SUBROUTINE  BCL4  (PTOTIN,RTOTIN,PREF,RREF) 

INTEGER  QOUT 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203),S(203) ,X(203) 
DIMENSION  XARRAYUoO) 

COMMON/ SUBS/P,R,U, A, S,X 

COMMON/ FUN 1 / G , PA , RA , UA , RB , RMU 

COMMON  XARRAY.Nl 

N1=N1+1 

QOUT=Nl/5 

PTOT=PTOTIN/PREF 

RTOT=RTOTIN/ RREF 

ATOT=  SQRT ( G*PTOT / RTOT ) 

STOT= ALOG ( PTOT/ RTOT**G ) 

U(1)=U(3) 

A(l)=SQRT(ATOT**2.-(G-l. )/2.*ABS(U(l))**2. ) 
AMACH=U(1)/A(1) 

IF(AMACH.LT.O.O)  GO  TO  60 

P(l)=PTOT/(l.+(G-l. )/2.*AMACH**2. )**(G/(G-1. )) 
R(1)=RT0T/U.  +  (G-1-  )/2-*AMACH**2.  )) 

S(l)=ALOG(P(l)/R(l)**G) 

GO  TO  50 

60  P(l)=PTOT/(l.+(G-l. )/2.*ABS(AMACH)**2. )**(G/(G-1. )) 
R( 1) =RTOT/ ( 1 . + (G- 1 . ) / 2 . *ABS (AMACH)**2 . )**( 1 • / (G- 1 . ) ) 
S(l)=ALOG(P(l)/R(l)**G) 

50  RETURN 
END 

SUBROUTINE  BCR4  (PTOTIN,RTOTIN,PREF,RREF) 

INTEGER  QOUT 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203) ,S(203)  ,X(203) 
DIMENSION  XARRAY(IOO) 

COMMON/ SUBS/P, R,U, A, S,X 

COMMON / FUN 1 / G , PA , RA , UA , RB , RMU 

COMMON  XARRAY.Nl 

N1=N1+1 

QOUT=Nl/25 

PTOT=PTOTIN/PREF 

RTOT=RTOTIN/RREF 

ATOT = SQRT ( G*PTOT / RTOT ) 

STOT= ALOG ( PTOT/ RTOT**G ) 

U(203)=U(201) 

A(203)=SQRT(ATOT**2. - (G- 1 . ) /2 . *ABS (U(203 ) )**2 . ) 
AMACH=U(203)/A(203) 
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IF(AMACH.LT.O.O)  GO  TO  60 

P(203)=PTOT/(1.+(G-1. )/2.*AMACH**2. )**(G/(G-1. )) 
R(203)=RTOT/(1.  +  (G-1.  )/2.*AMACH**2.  )**U./(G-1.  )) 
S(203)=ALOG(P(203)/R(203)**G) 

GO  TO  50 

60  P(203)=PTOT/(1.+(G-1. )/2.*ABS(AMACH)**2. )**(G/(G-1. )) 
rU03)=RTOT/U-  +  (G-1.  )/2.*ABS(AMACH)**2.  )**(!. /(G-1.  )) 
S(203)=ALOG(P(203)/R(203)**G) 

50  RETURN 
END 

SUBROUTINE  BCL5 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203) ,S(203) ,X(203) 

COMMON/ SUBS/P,R,U, A, S,X 

P(1)=P(3) 

R(1)=R(3) 

U(1)=U(3) 

A(1)=A(3) 

RETURN 

END 

SUBROUTINE  BCR5 

DIMENSION  P(203) ,R(203),U(203) ,A(203),S(203) ,X(203) 

COMMON/ SUBS/P,R,U, A, S,X 

P(203)=P(201) 

R(203)=R(201) 

U(203)=U(201) 

A(203)=A(201) 

RETURN 

END 

FUNCTION  WDP(II) 

DIMENSION  WNORM(12) ,IDIGT(12) 

COMMON/ SAMPLE/ WNORM , IDIGT 
IF  (lI.EQ.O)  GO  TO  10 
Ll=2 
L2=l 

DO  20  JJ=1,12 
IDIGT(JJ)=0 

WNORM ( J J ) = 1 . / FLOAT (LI** JJ ) 

20  CONTINUE 
WDP=0. 

RETURN 

10  DO  40  JJ=1,12 
Ll=2 
L2  =  l 

KJO=IDIGT(JJ) 

KJN=MOD((KJO+1) ,L1) 

IDIGT (JJ)=KJN 
IF  (KJO.LT.KJN)  GO  TO  50 
40  CONTINUE 
50  SUM=0. 

DO  60  JJ=1,12 

KNEW=MOD ( IDIGT ( J J )*L2 , LI ) 

SUM=  SUM+ FLO AT ( KNEW ) *WNORM ( J J ) 

60  CONTINUE 
WDP=SUM 
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RCM05180 

RCM05190 

RCM05200 

RCM05210 

RCM05220 

RCM05230 

RCM05240 

RCM05250 

RCM05260 

RCM05270 

RCM05280 

RCM05290 

RCM05300 

RCM05310 

RCM05320 

RCM05330 

RCM05340 

RCM05350 

RCM05360 

RCM05370 

RCM05380 

RCM05390 

RCM05400 

RCM05410 

RCM05420 


■ 


RETURN 

END 

SUBROUTINE  PLOTl(K) 

DIMENSION  X0RG(4) ,Y0RG(4) ,YMAX(4) ,YMIN(4) 

DATA  XORG/0.5,4.75,0.5,4.75/ 

DATA  YORG/0.5,0.5,4.75,4.75/ 

DATA  YMAX/3.50,3.0,0.5,2.0/ 

DATA  YMIN/0.5,0.0,-0.5,-2.0/ 

DO  10  1=1,4 

CALL  PHYSOR(XORG(I) ,YORG(I)) 

CALL  AREA2D(3.5,3.5) 

CALL  FRAME 

CALL  GRAF(0. , ’SCALE' , 1 . 0 ,YMIN(I ),' SCALE' ,YMAX(I)) 
CALL  ENDGR(O) 

10  CONTINUE 

CALL  PHYSOR(8.5,0.5) 

CALL  AREA2D(2.25,7.75) 

CALL  FRAME 

CALL  GRAF(0. , 'SCALE' ,1. ,0, 'SCALE' ,K) 

CALL  ENDGR(O) 

RETURN 

END 


RCM05430 

RCM05440 

RCM05450 

RCM05460 

RCM05470 

RCM05480 

RCM05490 

RCM05500 

RCM05510 

RCM05520 

RCM05530 

RCM05540 

RCM05550 

RCM05560 

RCM05570 

RCM05580 

RCM05590 

RCM05600 

RCM05610 

RCM05620 

RCM05630 

RCM05640 


SUBROUTINE  PLOT2(N,K)  RCM05650 

DIMENSION  X0RG(4),YORG(4),YMAX(4),YMIN(4),KNT(4),IYNAM(10)  RCM05660 

DIMENSION  PARRAY ( 100 ) , RARRAY ( 100 ) ,UARRAY UOO ) , SARRAY ( 100 ) , XARRAY ( 1RCM05 6 70 
*00)  RCM05680 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203)  RCM05690 

COMMON/SUBS/P, R,U, A, S,X  RCM05700 

COMMON  XARRAY  RCM05710 


DATA  XORG/0.5,4.75,0.5,4.75/ 

DATA  YORG/0.5,0.5,4.75,4.75/ 

DATA  YMAX/3.50,3.0,0.5,2.0/ 

DATA  YMIN/0.5,0.0,-0.5,-2.0/ 

DATA  KNT/1,4,6,9/ 

DATA  lYNAM/ ' PRES ’ , ’ SURE ' , ' $  ’ , ’ DENS ’ , ' ITY$ ' , ’ VELO ' , ’ CITY ' , ' $ 

*, 'ENTR' , 'OPY$' / 

DO  200  1=1,100 
PARRAY(I)=P (1*2+1) 

RARRAY(I)=R(I*2+1) 

UARRAY(I)=U( 1*2+1) 

SARRAY(I)=S(I*2+1) 

200  CONTINUE 


RCM05720 

RCM05730 

RCM05740 

RCM05750 

RCM05760 

’RCM05770 

RCM05780 

RCM05790 

RCM05800 

RCM05810 

RCM05820 

RCM05830 

RCM05840 


DO  300  1=1,4 

CALL  PHYSOR(XORG(I) ,YORG(I) ) 

CALL  AREA2D(3.5,3.5) 

CALL  XNAME( 'X' ,1) 

CALL  YNAME ( I YNAM ( KNT ( I ) ) , 10  0 ) 

CALL  GRAF(0. , 'SCALE' , 1 . 0 ,YMIN(I ),’ SCALE' ,YMAX(I)) 
IF(I.EQ.l)  CALL  SETCLR( 'YELLOW ) 

IF ( I . EQ . 2 )  CALL  SETCLR ( ’ CYAN ' ) 

IF(I.EQ.3)  CALL  SETCLR( 'RED' ) 

IF(I.EQ.4)  CALL  SETCLR ( 'MAGENTA' ) 

IF(N.EQ.K)  CALL  SETCLR ( 'WHITE' ) 

IF(I.EQ.l)  CALL  CURVE  (XARRAY, PARRAY, 100,0) 


RCM05850 

RCM05860 

RCM05870 

RCM05880 

RCM05890 

RCM05900 

RCM05910 

RCM05920 

RCM05930 

RCM05940 

RCM05950 

RCM05960 
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IF(I.EQ.2)  CALL  CURVE  (XARRAY.RARRAY, 100 , 0) 

RCM05970 

\  A 

IF(I.EQ.3)  CALL  CURVE  (XARRAY.UARRAY, 100,0) 

RCM05980 

IF(I.EQ.4)  CALL  CURVE  (XARRAY,SARRAY, 100 ,0) 

RCM05990 

CALL  ENDGR(O) 

RCM06000 

300 

CONTINUE 

RCM06010 

9m 

RETURN 

RCM06020 

m 

END 

RCM06030 

SUBROUTINE  GE ( SWL , SWR , N , TTOTAL , TIME , UEXMAX , PTOTIN , PREF ) 

RCM06040 

■  j 

INTEGER  SWL, SWR 

RCM06050 

si 

DIMENSION  P(203) ,R(203) ,U(203) ,A(203),S(203),X(203) 

RCM06060 

I 

COMMON/ SUBS/P,R,U, A, S,X 

RCM06070 

1  ^ 

C**CALCULATION  STARTS  AT  EXHAUST  PORT  OPENING.  SUBROUTINE  STRUCTURED 

RCM06080 

t  ^ 

C**ACCORDINGLY. 

RCM06090 

•  *  -  ■  / 

L 

IF((SWL.EQ.1).AND. (SWR.EQ.2))  GO  TO  10 

RCM06100 

.  *■*  i 

1  . 

IF((SWL.EQ.4).AND. (SWR.EQ.2))  GO  TO  30 

RCM06110 

-■  I'' ■ 

IF((SWL.EQ.4).AND. (SWR.EQ.l))  GO  TO  50 

RCM06120 

■;  t: 

IF( (SWL.EQ. 1) .AND. (SWR.EQ.l))  RETURN 

RCM06130 

»  •  •  • 

10 

PWALL=P(2) 

RCM06140 

“■  V, 

IF ( P WALL. LE. (PTOTIN /PREF))  GO  TO  20 

RCM06150 

jj  ^ 

RETURN 

RCM06160 

20 

SWL=4 

RCM06170 

F 

WRITE(6,74) 

RCM06180 

■b- 

WRITE (6, 75)  N, TTOTAL, TIME 

RCM06190 

■  .;-r. 

t 

RETURN 

RCM06200 

K 

30 

UEXIT=U(202) 

RCM06210 

li  k 

IF (UEXMAX. LT.UEXIT)  UEXMAX=UEXIT 

RCM06220 

•* 

1 1 

IF (UEXIT.lt. UEXMAX/ 2. )  GO  TO  40 

RCM06230 

RETURN 

RCM06240 

1  ■ 

40 

SWR=1 

RCM06250 

•'  .v. 

WRITE(6,76) 

RCM06260 

‘.■■V 

t 

WRITE (6, 75)  N, TTOTAL, TIME 

RCM06270 

t' 

RETURN 

RCM06280 

1  p 

50 

P1SH0K=P(2) 

RCM06290 

F  ? 

IF(P1SH0K.GT.PT0TIN/PREF)  GO  TO  60 

RCM06300 

RETURN 

RCM06310 

:  r 

60 

SWL=1 

RCM06320 

H 

WRITE(6,77) 

RCM06330 

1 

74 

WRITE (6, 75)  N, TTOTAL, TIME 

FORMAT (5X, 'INLET  PORT  OPENS  AT:’) 

RCM06340 

RCM06350 

*« 

k: — f 

75 

FORMAT(5X,I4,5X,2F14.7) 

RCM06360 

, 

76 

FORMAT ( 5X ,’ EXHAUST  PORT  CLOSES  AT:') 

RCM06370 

77 

FORMAT ( 5 X, 'INLET  PORT  CLOSES  AT:') 

RCM06380 

RETURN 

RCM06390 

END 

RCM06400 

1 

SUBROUTINE  SPCTRA(N , SWL , SWR , TIME , UEXMAX , PSEXIT , PSOUTl , PSOUT2 , PSOUTRCM06410 

*3 , JCOUNT , QPRINT , TTOTAL , KCOUNT ) 

RCM06420 

as 

INTEGER  SWL, SWR, QPRINT 

RCM06430 

.0.1 

t ' ' 

DIMENSION  P(203),R(203),U(203),A(203),S(203),X(203) 

RCM06440 

'S'  ■ 

fc 

COMMON/SUBS/P,R,U,A,S,X 

RCM06450 

;ov 

C*CALCULATION  STARTS  AT  HP  GAS  IN  PORT.  JCOUNT  IS  NUMBERED  ACCORDINGLY* 

*RCM06460 

•  '  < 

IF((SWL.EQ.1).AND. (SWR.EQ.3))  GO  TO  10 

RCM06470 

r 

IF((SWL.EQ.2) .AND. (SWR.EQ.3))  GO  TO  20 

RCM06480 

Kj 

IF((SWL.EQ. 2) .AND. (SWR.EQ.l))  GO  TO  30 

RCM06490 

IF((SWL.EQ. 5) .AND. (SWR.EQ.l))  GO  TO  40 

RCM06500 

H 

»,* 

iip^i 

IF((SWL.EQ.2) .AND. (SVm.EQ.5))  GO  TO  50  RCM06510 

IF((SWL.EQ.2).AND. (SWR.EQ.4))  GO  TO  60  RCM06520 

IF((SWL.EQ.1).AND. (SWR.EQ.4))  GO  TO  70  RCM06530 

10  IF(U(3) .LT.0.0)  GO  TO  11  RCM06540 

RETURN  RCM06550 

11  JCOUNT=JCOUNT+l  RCM06560 

PSEXIT=PS0UT1  RCM06570 

SWL=2  RCM06580 

WRITE(6,12)  RCM06590 

WRITE ( 6 , 1 3 ) N , TIME , S WL , S WR , J  COUNT  RCMO  6600 

12  FORMAT(5X, ’H.P.  AIR  OUT  PORT  OPENS  AT’)  RCM06610 

13  F0RMAT(5X,I4,5X.F9.7,5X,3I3)  RCM06620 

RETURN  RCM06630 

20  DO  26  1=5,199,2  RCM06640 

IF((R(I)-R(I+2)).GT.0.1)  GO  TO  22  RCM06650 

GO  TO  26  RCM06660 

22  XCNTCT=X(I)  RCM06670 

UCNTCT=U(I)  RCM06680 

TCNTCT=XCNTCT/ABS(UCNTCT)  RCM06690 

AHEAD=A(199)  RCM06700 

THEAD=1.0/A(199)  RCM06710 

IF(TCNTCT.LE.THEAD)  GO  TO  23  RCM06720 

RETURN  RCM06730 

23  JCOUNT=JCOUNT+l  RCM06740 

SWR=1  RCM06750 

WRITE(6,24)  RCM06760 

WRITE ( 6 , 25 )N , TIME , SWL , SWR , JCOUNT  RCM06  7  70 

24  FORMAT(5X, 'H.P.  GAS  IN  PORT  CLOSES  AT’)  RCM06780 

25  F0RMAT(5X,I4,5X,F9.7,5X,3I3)  RCM06790 

RETURN  RCM06800 

26  CONTINUE  RCM06810 

30  IF ( JCOUNT. EQ. 4)  GO  TO  80  RCM06820 

IF ( JCOUNT. EQ. 6)  GO  TO  90  RCM06830 

IF ( JCOUNT. EQ. 8)  GO  TO  100  RCM06840 

IF((R(2)-R(4)).GT.0.1)  GO  TO  31  RCM06850 

RETURN  RCM06860 

31  JCOUNT=JCOUNT+l  RCM06870 

SWL=5  RCM06880 

WRITE (6, 32)  RCM06890 

WRITE ( 6 , 3  3 ) N , TIME , S WL , SWR , JCOUNT  RCMO  6900 

32  FORMAT ( 5 X, ’HP  AIR  OUT  PORT  CLOSES  AND  TUNING  PORT  LI  OPENS  AT’)  RCM06910 

33  F0RMAT(5X,I4,5X,F9.7,5X,3I3)  RCM06920 

RETURN  RCM06930 

40  IF ( JCOUNT. EQ. 7)  GO  TO  110  RCM06940 

IF(U(3).GE.0.0)  GO  TO  41  RCM06950 

RETURN  RCM06960 

41  JCOUNT=JCOUNT+l  RCM06970 

PSEXIT=PS0UT2  RCM06980 

SWL=2  RCM06990 

WRITE (6, 42)  RCM07000 

WRITE ( 6 , 4  3 ) N , TIME , SWL , SWR , JCOUNT  RCMO  7010 

42  FORMAT (5X, ’TUNING  PORT  LI  CLOSES  AND  EXHAUST  PORT  El  OPENS  AT')  RCM07020 

43  F0RMAT(5X,I4,5X,F9.7,5X,3I3)  RCM07030 

RETURN  RCMO 7 040 
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80  IF(ABS(U(201)) .GT. .0001)  GO  TO  81 
RETURN 

81  JCOUNT=JCOUNT+l 
WRITE(6,82) 

WRITE ( 6 , 8  3 )N , TIME , SWL , SWR , JCOUNT 

82  FORMAT (5X, 'TUNING  PORT  R1  OPENS  AT’) 

83  F0RMAT(5X,IA,5X,F9.7,5X,3I3) 

RETURN 

50  IF ( JCOUNT. EQ. 9)  GO  TO  120 
IF(ABS(U(201)-U(3)) .LE. 0.001)  GO  TO  51 
RETURN 

51  JCOUNT= JCOUNT*! 

SWR=1 

WRITE (6, 52) 

WRITE ( 6 , 5  3 )N , TIME , SWL , SWR , JCOUNT 

52  FORMAT (5X, 'TUNING  PORT  R1  CLOSES  AT') 

53  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

90  THEAD1=X(201)/(ABS(U(3))*A(3)) 

KCOUNT=KCOUNT* 1 
IF(KCOUNT.EQ. 1)  TTOTl-TTOTAL 
IF(TTOTAL.GE. (TT0T1*THEAD1))  GO  TO  91 
RETURN 

91  JCOUNT= JCOUNT*! 

SWL=5 

WRITE(6,92) 

WRITE ( 6 , 93 )N , TIME , SWL , SWR , JCOUNT 

92  FORMAT (5X, 'EXHAUST  PORT  El  CLOSES  AND  TUNING  PORT  L2  OPENS  AT’) 

93  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

110  IF(U(3).GE.0.0)  GO  TO  111 
RETURN 

111  JCOUNT= JCOUNT*! 

PSEXIT=PSOUT3 

CUT  =9 

WRITE (6, 112) 

WRITE ( 6 , 113 )N , TIME , SWL , SWR , JCOUNT 

112  FORMAT (5X, ’TUNING  PORT  L2  CLOSES  AND  EXHAUST  PORT  E2  OPENS  AT') 

113  F0RMAT(5X,IA,5X,F9.7,5X,3I3) 

RETURN 

100  IF(ABS(U(201)) .GT. .0001)  GO  TO  101 
RETURN 

101  J COUNT = JCOUNT*! 

SWR=  5 

WRITE(6,102) 

WRITE ( 6 , 103 )N , TIME , SWL , SWR , JCOUNT 

102  FORMAT ( 5 X, 'TUNING  PORT  R2  OPENS  AT’) 

103  F0RMAT(5X,I4,5X,F9.7,5X,3I3) 

RETURN 

120  IF(ABS(U(201)-U(3)) .LE. 0.0001)  GO  TO  121 
RETURN 

121  JCOUNT= JCOUNT*! 

SWR=4 


RCM07050 

RCM07060 

RCM07070 

RCM07080 

RCM07090 

RCM07100 

RCM07110 

RCM07120 

RCM07130 

RCM07140 

RCM07150 

RCM07160 

RCM07170 

RCM07180 

RCM07190 

RCM07200 

RCM07210 

RCM07220 

RCM07230 

RCM07240 

RCM07250 

RCM07260 

RCM07270 

RCM07280 

RCM07290 

RCM07300 

RCM07310 

RCM07320 

RCM07330 

RCM07340 

RCM07350 

RCM07360 

RCM07370 

RCM07380 

RCM07390 

RCM07400 

RCM07410 

RCM07420 

RCM07430 

RCM07440 

RCM07450 

RCM07460 

RCM07470 

RCM07480 

RCM07490 

RCM07500 

RCM07510 

RCM07520 

RCM07530 

RCM07540 

RCM07550 

RCM07560 

RCM07570 

RCM07580 


WRITE(6,122)  RCM07590 

VrRITE(6,l23)N,TIME,SWL,SWR,  JCOUNT  RCMO  7  6  0 0 

122  FORMAT (5X, 'TUNING  PORT  R2  CLOSES  AND  L.P.  AIR  INLET  OPENS  AT’)  RCM07610 

123  FORMAT(5X,I4,5X,F9.7,5X,3I3)  RCM07620 

RETURN  RCMO 7 6 30 

60  IF((R(4)-R(2)).GT.0.1)  GO  TO  61  RCM07640 

RETURN  RCM07650 

61  JCOUNT=JCOUNT+l  RCM07660 

SWL=1  RCMO 7 6 70 

WRITE(6,62)  RCM07680 

WRITE ( 6 , 6  3 ) N , TIME , SWL , S WR , JCOUNT  RCMO  7690 

62  FORMAT (5X, 'EXHAUST  PORT  E2  CLOSES  AT')  RCM07700 

63  FORMAT(5X,I4,5X,F9.7,5X,3I3)  RCM07710 

RETURN  RCMO  7  7  2  0 

70  IF(U(201) .GE.0.0)  GO  TO  71  RCM07730 

RETURN  RCMO 7 740 

71  JCOUNT=0  RCM07750 

SWR= 1  RCMO  7760 

WRITE(6,72)  RCM07770 

WRITE ( 6 , 7  3 ) N , TIME , S WL , S WR , JCOUNT  RCMO  7780 

72  FORMAT ( 5 X, 'CYCLE  COMPLETED.')  RCM07790 

73  FORMAT(5X,I4,5X,F9.7,5X,3I3)  RCM07800 

RETURN  RCMO  7810 

Fwn  RCMO 7 8 20 
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APPENDIX  B 
PROGRAM  RCM 

B. 1 ,  Program  Description 
B.l.l.  Computational  Grid 

The  computational  region  is  divided  into  \00  cells;  the  solution  grid 
points  are  odd  numbered,  e.g.,  3,  5,  7  201  with  1  and  203  being  the 

points  where  the  boundary  conditions  are  specified.  The  even  numbered  points, 
2,  4,  6  ...,  202  are  intermediate  locations  where  solutions  are  stored  before 
being  assigned  to  the  solution  grid  points.  See  Fig.  (2). 

B.1.2.  Data  Input 

Data  for  various  ports  (exhaust,  inlet,  etc.)  is  specified  in 
dimensional  form  in  S.I.  units  (Pascal  (N/m^)  for  pressure,  kg/m^  for 
density,  m/s  for  velocity  etc.).  Reference  values  are  also  specified  in 
like  manner.  See  lines  RCM00210  through  00250. 

Initial  data  Is  specified  through  a  call  to  an  appropriate 
subroutine,  depending  on  where  the  calculation  is  started  for  a  particular 
wave  diagram.  For  the  example  given  in  section  II  on  the  Spectra  Technology 
wave  diagram,  the  computation  is  started  at  the  point  when  the  high  pressure 
gas  inlet  port  just  opens.  The  call  for  initial  data  is  made  to  subroutine 
1NIT3R,  which  prescribes  data  consistent  with  a  solid  wall  boundary  at  the 
left  and  a  'piston'  inflow  boundary  at  the  right. 


L' 


B.I.3.  Non-dlmensionallzatlon 

Non-dlmenslonallzation  is  carried  out  in  lines  00540  through  00610 
with  entropy  defined  as 

S  -  in  (£y) 

P  ^ 


B.l 


u 


ii.V. 

m 


vv 

mt 


Note  that  velocities  are  all  referred  to  a  reference  sonic  velocity  defined 


by 


aref  ” 


Pref 

Pref 


B.1.4.  Structure 

The  main  program  loop  starts  at  line  00630,  for  the  number  of  time 
steps  specified.  The  time  step  is  computed  according  to  the  appropriate  CFL 
condition  for  the  method,  and  a  random  number  for  the  time  step  is  generated 
by  a  call  to  the  function  subroutine  WDP. 

A  secondary  loop  to  define  the  sequence  of  local  Rlemann  problems  for 
the  time  step  is  set  up  at  line  00750.  For  each  Rlemann  problem  defined,  a 
call  is  made  to  subroutine  GLIMM  which  1)  solves  the  Rlemann  problem,  and  ii) 
samples  the  solution  using  the  random  number  generated.  The  subroutine  then 
returns  the  sampled  solution  as  the  parameters  PGLIM,  RGLIM,  UGLIM  for  the 
pressure,  density  and  velocity  respectively.  These  solutions  are  initially 
stored  in  the  even  numbered  intermediate  locations  on  the  grid,  and  are  then 
assigned  to  either  the  left  or  the  right  solution  grid  point  depending  on 
whether  the  random  number  was  in  the  negative  or  the  positive  half  of  the 
interval  respectively. 

A  call  is  then  made  to  one  of  the  modular  subroutines  structured  for 
particular  types  of  wave  diagrams,  lines  01050-01080,  and  the  others  are 
commented  out. 

Boundary  conditions  are  invoked  after  the  call  to  the  modular 
subroutines  which  return  the  proper  values  of  the  switches  SWL  and  SWR.  The 
structure  of  the  boundary  condition  subroutines  is  described  in  section  II. 
This  sequence  completes  one  pass  through  the  main  loop  and  the  process  is 
repeated  for  the  number  of  time  steps  specified. 


B.2.  Example  Use  of  Program  RCM 
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The  program  Is  set  up  In  the  following  steps: 

i)  Line  00150  -  output  device  designation.  See  B.3. 

11)  Line  00190  -  specify  the  number  of  time  steps,  k  ,  and  the 

switches  SWL  and  SUR  consistent  with  where  the 
computation  Is  to  be  started. 

iii)  Lines  00210  -  prescribe  flow  data  for  various  ports  in 
through  00250  dimensional  form.  See  list  of  variables  for 

explanation  of  variable  names. 

iv)  Lines  00490  -  invoke  the  proper  initial  data  subroutine  and 
through  00530  comment  out  the  rest.  See  list  of  subroutines 

for  explanation  of  subroutine,  function  subprogram 
names. 

v)  Line  00660  -  set  the  interval  for  number  of  time  steps  at  which 

a  plot  of  the  flow  parameters  is  required. 

vi)  Lines  01050  -  user  supplied  modular  subroutine  for  particular 
through  01080  wave  diagram  to  be  computed.  Comment  out  the  rest 

vii)  Line  01190  -  call  to  plotting  routine  should  be  consistent  with 

interval  specified  in  line  0660. 

viii)  Lines  02650  -  identify  proper  subroutine  to  prescribe  initial 
through  03700  data  (consistent  with  iv),  and  specify  the  data  in 

the  subroutine  in  dimensional  form. 

ix)  Lines  05470  -  specify  plotting  parameters,  viz.,  origins  of 
through  05990  plots,  scales,  number  of  points  to  be  plotted, 

color  of  plots,  etc.  Facility  dependent. 

The  subroutines  PLOTl  and  PL0T2  given  in  the 
listing  are  structured  for  DISSPLA  software 
installed  in  the  facility  at  NPGS. 

x)  Lines  06040  -  user  supplied  modular  subroutine  for  wave  diagram 
through  07820  to  be  computed. 


i.. 

B.3.  Execution 

.y  The  program  is  run  in  an  interactive  mode  and  is  invoked  through  a  call 

to  DISSPLA,  available  on  most  mainframes.  After  compiling  the  program 


B.3 


(FORTRAN  H  Extended  compiler),  the  following  command  executes  it: 


DISSPLA  filename 

If  working  at  stations  equipped  with  dual  screens,  the  command  on  line  150  can 
be  of  the  type 

CALL  TEK618  +  Tektronix  screen 

If  working  on  a  non-graphics  terminal,  or  a  single  screen  station,  this 
should  be  changed  to 

CALL  COMPRS 

which  generates  a  'DISSPLA  METAFILE'  to  be  routed  later  to  either  a  screen  or 
a  plotter,  e.g.,  VRSTEC,  IBM79,  TEK618,  etc.  Once  generated,  the  metafile  can 
be  accessed  and  routed  by  the  command 

DISSPOP  device  designation 

These  are  facility  dependent  commands  and  should  be  modified  accordingly. 


B.4.  List  of  Important  Variables  (In  Alphabetical  Order) 


A 

AHEAD 

AMACH 

AL 

AR 

AREF 

ASTAR 

CFLNUM 

DT 

DX 

EPS 

G 

IDIGT 

II 

JCOUNT 

K 

KCOUNT 

N 

N1 


sonic  velocity 

sonic  speed  of  head  wave  of  rarefaction  fan 
Mach  number 

left  side  sonic  speed  value  for  RP 
right  side  sonic  speed  value  for  RP 
reference  speed  of  sound 

speed  of  sound  in  'starred'  state  of  RP  solution 
(see  Fig.  3) 

CFL  number  for  time  step  determination 

time  step 

grid  cell  width 

small  number  for  pressure  iteration  in  RP  solver 
ratio  of  specific  heats,  y 
see  WNORM 

argument  used  in  function  subprogram  PHI  equal  to  either 
0  or  1 
counter 

number  of  time  steps 
counter 

counter  for  time  steps 
counter 


p 

PA 

PGLIM 

PL 

PR 

PREF 

PSEXIT 

PSINL 

PSINR 

PSOUTn 

PSTAR 

PTOTIN 

QPRINT 

QSTOP 

R 

RA,RB 

RGLIM 

RINL 

RINR 

RL 

RMU 

RR 

RREF 

RTOTIN 

S 

SWL 

SWR 

TCNTCT 

THEAD 

TIME 
TIME REF 
TTOTAL 
U 

UA 

UCNTCT 

UEXMAX 

UGLIM 

UL 

UR 

USTAR 

WDP 

WL 

WR 

WNORM 

X 

XCNTCT 
XI, XII 


pressure 

flow  parameter  describing  'a'  state  in  transition  functions 

pressure  value  returned  by  subroutine  GLIMM 

left  side  pressure  value  for  RP 

right  side  pressure  value  for  RP 

reference  pressure 

static  pressure  at  exit  or  outlet  port 
static  pressure  for  incoming  'piston'  flow  on  left  side 
static  pressure  for  incoming  'piston'  flow  on  right  side 
n  =  1,2,3  -  exit  static  pressures  for  cycles  with  more  than 
one  exhaust  port 

pressure  in  'starred'  state  of  RP  solution  (see  Fig.  3) 
total  pressure  for  isentropic  Inflow 
specification  of  interval  size  for  output 
maximum  number  of  iterations  for  solution  of  Rlemann 
problem,  (RP) 
density 

flow  parameters  describing  'a'  and  'b'  states  in  transition 
functions 

density  returned  by  subroutine  GLIMM 

static  density  for  incoming  'piston'  flow  on  left  side 

static  density  for  incoming  'piston'  flow  on  right  side 

left  side  density  for  RP 

function  of  y 

right  side  density  for  RP 

reference  density 

total  density  for  Isentropic  inflow 
entropy 

switch  for  left  boundary 
switch  for  right  boundary 

time  taken  by  contact  surface  to  travel  a  certain  distance 
time  taken  by  head  wave  of  expansion  to  travel  a  certain 
distance 

real  time  in  seconds 
reference  time 

cumulative  non-dimensional  time  for  number  of  time  steps 
velocity 

flow  parameter  for  'a'  state  in  transition  functions 
velocity  of  contact  surface 

maximum  velocity  occurring  at  an  outflow  boundary 
velocity  returned  by  subroutine  GLIMM 
left  side  velocity  for  RP 
right  side  velocity  for  RP 

velocity  in  'starred'  state  of  RP  solution  (see  Fig,  3) 
value  returned  by  random  number  generator  subprogram 
left  shock  wave  velocity 
right  shock  wave  velocity 

variable  used  in  random  number  generator  subprogram 

space  dimension 

location  of  contact  surface 

random  numbers  scaled  to  grid  cell 
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XREF 

Y 

Z 

ZETA 


reference  length 

argument  used  in  function  subprogram  PHI  equal  to  PSTAR 
argument  used  in  function  subprogram  PHI  equal  to  either 
PL  or  PR 

dummy  variable  (for  initialization  purposes  in  random 
number  generator ) 


B. 5.  List  of  Subroutines,  Function  Subprograms 


B.5.1.  Subroutines 


INITl 

-  prescribes  initial  data  corresponding 
e.g.,  shock-tube  problem 

to 

SWL=1, 

SWR=1; 

INIT2L 

-  prescribes 

initial 

data 

corresponding 

to 

SWL=2, 

SWR=1 

INIT2R 

-  prescribes 

initial 

data 

corresponding 

to 

SWL=1, 

SWR=2 

INIT3L 

-  prescribes 

Initial 

data 

corresponding 

to 

SWL=3, 

SWR=1 

INIT3R 

-  prescribes 

initial 

data 

corresponding 

to 

SWL=1, 

SWR=3 

PL0T1,2  -  graphics  subroutines 
GLIMM 


GE 


solves  the  Riemann  problem,  samples  the  solution  and 
returns  values  for  flow  parameters 

modular  user  supplied  subroutine  to  simulate  wave  diagram 
of  General  Electric  Wave  Engine 
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DETON  -  modular  user  supplied  subroutine  to  simulate  evacuation  of 
detonation  chamber 

SPCTRA  -  modular  user  supplied  subroutine  to  simulate  wave  diagram 
of  Spectra  Technology's  Pressure  Exchanger 

BCLl  -  prescribes  boundary  conditions  (BC's)  corresponding  to 
SWL*1,  i.e.,  solid  wall  on  left  side 

BCL2  -  prescribes  BC's  corresponding  to  SWL=2,  i.e.,  outflow  at 
constant  static  pressure  on  left  side 

BCL3  -  prescribes  BC's  corresponding  to  SWL=3,  i.e.,  'piston' 
Inflow  on  left  side 
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BCL4  -  prescribes  BC's  corresponding  to  SWL“4,  i.e.,  Isentroplc 
inflow  from  reservoir  on  left  side 


BCL5  -  prescribes  BC's  corresponding  to  SWL=5,  i.e.,  wave  'tuning' 
on  left  side 
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BCRl ,  BCR2,  BCR3,  BCR4,  BCR5  -  prescribe  BC's  corresponding  to 
SWR*1 ,2,3,4,5  respectively  on  right  side 


Function  Subprograms 

PHI(y,z)  -  required  in  Iteration  procedure  for  solution  of  RP 

PHIl(PB)  -  describes  shock  transition  function,  (Pa^Pb)  » 

two  states  a  and  b  connected  by  a  shock  wave  (see 
Ref.  6,  Ch.  Ill) 

PSI(PB)  -  describes  rarefaction  transition  function,  'Pa(pb)  » 
for  two  states  a  and  b  connected  by  a  rarefaction  wave 
(see  Ref.  6,  Ch.  Ill) 


WDP(II)  -  generates  a  random  number  in  a  van  der  Corput  sequence 

each  time  it  is  invoked.  Note  that  it  needs  to  be  called 
once  from  outside  the  main  loop  by  specifying  an  argument 
11=1  to  initialize  IDIGT  and  WNORM,  returning  a  value 
of  0  for  the  dummy  variable  ZETA,  and  then  a  second  time 
from  within  the  main  loop  with  an  argument  11=0  to  return 
a  value  which  is  the  random  number. 


